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

jengelh-datetimespielwiese
Last change on this file since c81423 was c81423, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
better condition for running "m_checker.Initialize(m_syzLeads);" in SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms add: CReducerFinder::DebugPrint() for hash output
  • Property mode set to 100644
File size: 27.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  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__ && !__IGNORETAILS__ && (IDELEMS(m_syzLeads) > 0) && !((IDELEMS(m_syzLeads) == 1) && (m_syzLeads->m[0] == NULL)))
622  {
623    m_LS = m_syzLeads;
624    m_checker.Initialize(m_syzLeads);
625#ifndef NDEBUG   
626    if( __DEBUG__ )
627    {
628      const ring& r = m_rBaseRing;
629      PrintS("SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms: \n");
630      PrintS("m_syzLeads: \n");
631      dPrint(m_syzLeads, r, r, 1);
632      PrintS("m_checker.Initialize(m_syzLeads) => \n");
633      m_checker.DebugPrint();
634    }
635#endif
636    assume( m_checker.IsNonempty() ); // TODO: this always fails... BUG????
637  }
638}
639
640poly SchreyerSyzygyComputation::SchreyerSyzygyNF(poly syz_lead, poly syz_2) const
641{
642// return FROM_NAMESPACE(INTERNAL, _SchreyerSyzygyNF(syz_lead, syz_2, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
643// poly _SchreyerSyzygyNF(poly syz_lead, poly syz_2,
644//                       ideal L, ideal T, ideal LS,
645//                       const ring r,
646//                       const SchreyerSyzygyComputationFlags attributes)
647// {
648
649  assume( !__IGNORETAILS__ );
650
651  const ideal& L = m_idLeads;
652  const ideal& T = m_idTails;
653  const ring& r = m_rBaseRing;
654
655//   const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
656//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
657//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
658//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
659//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
660//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
661
662  assume( syz_lead != NULL );
663  assume( syz_2 != NULL );
664
665  assume( L != NULL );
666  assume( T != NULL );
667
668  assume( IDELEMS(L) == IDELEMS(T) );
669
670  int  c = p_GetComp(syz_lead, r) - 1;
671
672  assume( c >= 0 && c < IDELEMS(T) );
673
674  poly p = leadmonom(syz_lead, r); // :( 
675  poly spoly = pp_Mult_qq(p, T->m[c], r);
676  p_Delete(&p, r);
677
678
679  c = p_GetComp(syz_2, r) - 1;
680  assume( c >= 0 && c < IDELEMS(T) );
681
682  p = leadmonom(syz_2, r); // :(
683  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
684  p_Delete(&p, r);
685
686  poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
687
688  while (spoly != NULL)
689  {
690    poly t = m_div.FindReducer(spoly, NULL, m_checker);
691
692    p_LmDelete(&spoly, r);
693
694    if( t != NULL )
695    {
696      p = leadmonom(t, r); // :(
697      c = p_GetComp(t, r) - 1;
698
699      assume( c >= 0 && c < IDELEMS(T) );
700
701      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
702
703      p_Delete(&p, r);
704
705      tail = p_Add_q(tail, t, r);
706    }   
707  }
708
709  return tail;
710}
711
712
713poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
714{
715  assume( !__IGNORETAILS__ );
716
717  const ideal& L = m_idLeads;
718  const ideal& T = m_idTails;
719//  const ideal& LS = m_LS;
720  const ring& r = m_rBaseRing;
721//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
722
723//  return FROM_NAMESPACE(INTERNAL, _TraverseTail(multiplier, tail, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
724// poly _TraverseTail(poly multiplier, poly tail,
725//                    ideal L, ideal T, ideal LS,
726//                    const ring r,
727//                    const SchreyerSyzygyComputationFlags attributes)
728// {
729
730
731//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
732//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
733//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
734//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
735//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
736
737  assume( multiplier != NULL );
738
739  assume( L != NULL );
740  assume( T != NULL );
741
742  poly s = NULL;
743
744  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
745    for(poly p = tail; p != NULL; p = pNext(p))   // iterate over the tail
746      s = p_Add_q(s, ReduceTerm(multiplier, p, NULL), r); 
747
748  return s;
749}
750
751
752
753
754poly SchreyerSyzygyComputation::ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck) const
755{
756  assume( !__IGNORETAILS__ );
757
758  const ideal& L = m_idLeads;
759  const ideal& T = m_idTails;
760//  const ideal& LS = m_LS;
761  const ring& r = m_rBaseRing;
762//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
763
764//  return FROM_NAMESPACE(INTERNAL, _ReduceTerm(multiplier, term4reduction, syztermCheck, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
765// poly _ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
766//                 ideal L, ideal T, ideal LS,
767//                 const ring r,
768//                 const SchreyerSyzygyComputationFlags attributes)
769
770
771
772//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
773//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
774//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
775//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
776//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
777
778  assume( multiplier != NULL );
779  assume( term4reduction != NULL );
780
781
782  assume( L != NULL );
783  assume( T != NULL );
784
785  // assume(r == currRing); // ?
786
787  // simple implementation with FindReducer:
788  poly s = NULL;
789
790  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
791  {
792    // NOTE: only LT(term4reduction) should be used in the following:
793    poly product = pp_Mult_mm(multiplier, term4reduction, r);
794    s = m_div.FindReducer(product, syztermCheck, m_checker);
795    p_Delete(&product, r);
796  }
797
798  if( s == NULL ) // No Reducer?
799    return s;
800
801  poly b = leadmonom(s, r);
802
803  const int c = p_GetComp(s, r) - 1;
804  assume( c >= 0 && c < IDELEMS(T) );
805
806  const poly tail = T->m[c];
807
808  if( tail != NULL )
809    s = p_Add_q(s, TraverseTail(b, tail), r); 
810
811  return s;
812}
813
814
815
816
817
818BEGIN_NAMESPACE_NONAME
819   
820static inline int atGetInt(idhdl rootRingHdl, const char* attribute, long def)
821{
822  return ((int)(long)(atGet(rootRingHdl, attribute, INT_CMD, (void*)def)));
823}
824
825END_NAMESPACE   
826
827SchreyerSyzygyComputationFlags::SchreyerSyzygyComputationFlags(idhdl rootRingHdl):
828#ifndef NDEBUG
829    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", TRUE) ),
830#else
831    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", FALSE) ),
832#endif
833//    __SYZCHECK__( (BOOLEAN)atGetInt(rootRingHdl, "SYZCHECK", __DEBUG__) ),
834    __LEAD2SYZ__( (BOOLEAN)atGetInt(rootRingHdl, "LEAD2SYZ", 1) ), 
835    __TAILREDSYZ__( (BOOLEAN)atGetInt(rootRingHdl, "TAILREDSYZ", 1) ), 
836    __HYBRIDNF__( (BOOLEAN)atGetInt(rootRingHdl, "HYBRIDNF", 0) ),
837    __IGNORETAILS__( (BOOLEAN)atGetInt(rootRingHdl, "IGNORETAILS", 0) ),
838    m_rBaseRing( rootRingHdl->data.uring )
839{   
840  if( __DEBUG__ )
841  {
842    PrintS("SchreyerSyzygyComputationFlags: \n");
843    Print("        DEBUG: \t%d\n", __DEBUG__);
844//    Print("   SYZCHECK  : \t%d\n", __SYZCHECK__);
845    Print("     LEAD2SYZ: \t%d\n", __LEAD2SYZ__);
846    Print("   TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
847    Print("  IGNORETAILS: \t%d\n", __IGNORETAILS__);
848   
849  }
850
851  // TODO: just current setting!
852  assume( rootRingHdl == currRingHdl );
853  assume( rootRingHdl->typ == RING_CMD );
854  assume( m_rBaseRing == currRing );
855  // move the global ring here inside???
856}
857
858   
859
860CReducerFinder::CLeadingTerm::CLeadingTerm(unsigned int _label,  const poly _lt, const ring R):
861    m_sev( p_GetShortExpVector(_lt, R) ),  m_label( _label ),  m_lt( _lt )
862{ }
863
864
865CReducerFinder::~CReducerFinder()
866{
867  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++ )
868  {
869    const TReducers& v = it->second;
870    for(TReducers::const_iterator vit = v.begin(); vit != v.end(); vit++ )
871      delete const_cast<CLeadingTerm*>(*vit);
872  }
873}
874
875
876void CReducerFinder::Initialize(const ideal L)
877{
878  assume( m_L == NULL || m_L == L );
879  if( m_L == NULL )
880    m_L = L;
881
882  assume( m_L == L );
883 
884  if( L != NULL )
885  {
886    const ring& R = m_rBaseRing;
887    assume( R != NULL );
888   
889    for( int k = IDELEMS(L) - 1; k >= 0; k-- )
890    {
891      const poly a = L->m[k]; // assume( a != NULL );
892
893      // NOTE: label is k \in 0 ... |L|-1!!!
894      if( a != NULL )
895        m_hash[p_GetComp(a, R)].push_back( new CLeadingTerm(k, a, R) );
896    }
897  }
898}
899
900CReducerFinder::CReducerFinder(const ideal L, const SchreyerSyzygyComputationFlags& flags):
901    SchreyerSyzygyComputationFlags(flags),
902    m_L(const_cast<ideal>(L)), // for debug anyway
903    m_hash()
904{
905  assume( flags.m_rBaseRing == m_rBaseRing );
906  if( L != NULL )
907    Initialize(L);
908}
909
910
911bool CReducerFinder::IsDivisible(const poly product) const
912{
913  const ring& r = m_rBaseRing;
914 
915  const long comp = p_GetComp(product, r);
916  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
917
918  assume( comp >= 0 );
919
920  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
921
922  if( it == m_hash.end() )
923    return false;
924
925  assume( m_L != NULL ); 
926
927  const TReducers& reducers = it->second;
928
929  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
930  {
931    const poly p = (*vit)->m_lt;
932
933    assume( p_GetComp(p, r) == comp );
934
935    const int k = (*vit)->m_label;
936
937    assume( m_L->m[k] == p );
938
939    const unsigned long p_sev = (*vit)->m_sev;
940
941    assume( p_sev == p_GetShortExpVector(p, r) );     
942
943    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
944      continue;
945
946    if( __DEBUG__ )
947    {
948      Print("_FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", k+1);
949      dPrint(p, r, r, 1);
950    }
951
952    return true;
953  }
954
955  return false;
956}
957
958
959#ifndef NDEBUG
960void CReducerFinder::DebugPrint() const
961{
962  const ring& r = m_rBaseRing;
963
964  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++)
965  {
966    Print("Hash Key: %d, Values: \n", it->first);
967    const TReducers& reducers = it->second;
968
969    for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
970    {
971      const poly p = (*vit)->m_lt;
972
973      assume( p_GetComp(p, r) == it->first );
974
975      const int k = (*vit)->m_label;
976
977      assume( m_L->m[k] == p );
978
979      const unsigned long p_sev = (*vit)->m_sev;
980
981      assume( p_sev == p_GetShortExpVector(p, r) );
982
983      Print("L[%d]: ", k); dPrint(p, r, r, 0); Print("SEV: %dl\n", p_sev);     
984    }
985  }
986}
987#endif
988
989
990poly CReducerFinder::FindReducer(const poly product, const poly syzterm, const CReducerFinder& syz_checker) const
991{
992  const ring& r = m_rBaseRing;
993
994  assume( product != NULL );
995
996  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
997
998  long c = 0;
999
1000  if (syzterm != NULL)
1001    c = p_GetComp(syzterm, r) - 1;
1002
1003  assume( c >= 0 && c < IDELEMS(L) );
1004
1005  if (__DEBUG__ && (syzterm != NULL))
1006  {
1007    const poly m = L->m[c];
1008
1009    assume( m != NULL ); assume( pNext(m) == NULL );
1010
1011    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
1012    assume( p_EqualPolys(lm, product, r) );
1013
1014    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
1015    //  def @@product = leadmonomial(syzterm) * L[@@r];
1016
1017    p_Delete(&lm, r);   
1018  }
1019
1020  const long comp = p_GetComp(product, r);
1021  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
1022
1023  assume( comp >= 0 );
1024
1025//   for( int k = IDELEMS(L)-1; k>= 0; k-- )
1026//   {
1027//     const poly p = L->m[k];
1028//
1029//     if ( p_GetComp(p, r) != comp )
1030//       continue;
1031//
1032//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
1033 
1034   // looking for an appropriate diviser p = L[k]...
1035  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
1036
1037  if( it == m_hash.end() )
1038    return NULL;
1039
1040  assume( m_L != NULL );
1041
1042  const TReducers& reducers = it->second;
1043
1044  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
1045
1046  const poly q = p_New(r); pNext(q) = NULL;
1047
1048  if( __DEBUG__ )
1049    p_SetCoeff0(q, 0, r); // for printing q
1050 
1051  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1052  {
1053    const poly p = (*vit)->m_lt;
1054
1055    assume( p_GetComp(p, r) == comp );
1056
1057    const int k = (*vit)->m_label;
1058
1059    assume( L->m[k] == p );
1060
1061    const unsigned long p_sev = (*vit)->m_sev;
1062
1063    assume( p_sev == p_GetShortExpVector(p, r) );     
1064
1065    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
1066      continue;     
1067
1068//     // ... which divides the product, looking for the _1st_ appropriate one!
1069//     if( !p_LmDivisibleByNoComp(p, product, r) ) // included inside  p_LmShortDivisibleBy!
1070//       continue;
1071
1072    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
1073    p_SetComp(q, k + 1, r);
1074    p_Setm(q, r);
1075
1076    // cannot allow something like: a*gen(i) - a*gen(i)
1077    if (syzterm != NULL && (k == c))
1078      if (p_ExpVectorEqual(syzterm, q, r))
1079      {
1080        if( __DEBUG__ )
1081        {
1082          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1083          dPrint(syzterm, r, r, 1);
1084        }   
1085
1086        continue;
1087      }
1088
1089    // while the complement (the fraction) is not reducible by leading syzygies
1090    if( to_check && syz_checker.IsDivisible(q) ) 
1091    {
1092      if( __DEBUG__ )
1093      {
1094        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
1095      }
1096     
1097      continue;
1098    }
1099
1100    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
1101    return q;
1102  }
1103
1104  p_LmFree(q, r);
1105
1106  return NULL;
1107}
1108
1109
1110
1111CLCM::CLCM(const ideal& L, const SchreyerSyzygyComputationFlags& flags):
1112    SchreyerSyzygyComputationFlags(flags), std::vector<bool>(),
1113    m_compute(false), m_N(rVar(flags.m_rBaseRing))
1114{
1115  const ring& R = m_rBaseRing;
1116  assume( flags.m_rBaseRing == R );
1117  assume( R != NULL );
1118
1119  assume( L != NULL );
1120
1121  if( __TAILREDSYZ__ && !__HYBRIDNF__ && (L != NULL))
1122  {
1123    const int l = IDELEMS(L);
1124
1125    assume( l > 0 );
1126
1127    resize(l, false);
1128
1129    for( int k = l - 1; k >= 0; k-- )
1130    {
1131      const poly a = L->m[k]; assume( a != NULL );
1132
1133      for (unsigned int j = m_N; j > 0; j--)
1134        if ( !(*this)[j] )
1135          (*this)[j] = (p_GetExp(a, j, R) > 0);
1136    }
1137
1138    m_compute = true;   
1139  }
1140}
1141
1142
1143bool CLCM::Check(const poly m) const
1144{
1145  assume( m != NULL );
1146  if( m_compute && (m != NULL))
1147  { 
1148    const ring& R = m_rBaseRing;
1149
1150    assume( __TAILREDSYZ__ && !__HYBRIDNF__ );
1151
1152    for (unsigned int j = m_N; j > 0; j--)
1153      if ( (*this)[j] )
1154        if(p_GetExp(m, j, R) > 0)
1155          return true;
1156
1157    return false;
1158
1159  } else return true;
1160}
1161
1162
1163
1164
1165END_NAMESPACE               END_NAMESPACE_SINGULARXX
1166
1167
1168// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
Note: See TracBrowser for help on using the repository browser.