source: git/dyn_modules/syzextra/mod_main.cc @ f12f9a

spielwiese
Last change on this file since f12f9a was f12f9a, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
separation of low-level functions: ComputeLeadingSyzygyTerms & Compute2LeadingSyzygyTerms
  • Property mode set to 100644
File size: 55.3 KB
Line 
1
2
3
4
5#include <kernel/mod2.h>
6
7#include <omalloc/omalloc.h>
8
9#include <misc/intvec.h>
10#include <misc/options.h>
11
12#include <coeffs/coeffs.h>
13
14#include <polys/PolyEnumerator.h>
15
16#include <polys/monomials/p_polys.h>
17#include <polys/monomials/ring.h>
18// #include <kernel/longrat.h>
19#include <kernel/GBEngine/kstd1.h>
20
21#include <kernel/polys.h>
22
23#include <kernel/GBEngine/syz.h>
24
25#include <Singular/tok.h>
26#include <Singular/ipid.h>
27#include <Singular/lists.h>
28#include <Singular/attrib.h>
29
30#include <Singular/ipid.h> 
31#include <Singular/ipshell.h> // For iiAddCproc
32
33// extern coeffs coeffs_BIGINT
34
35#include "singularxx_defs.h"
36#include "DebugPrint.h"
37#include "myNF.h"
38
39
40#include <Singular/mod_lib.h>
41
42
43#if GOOGLE_PROFILE_ENABLED
44#include <google/profiler.h>
45#endif // #if GOOGLE_PROFILE_ENABLED
46
47
48#include <stdio.h>
49#include <stdlib.h>
50#include <string.h>
51
52
53#ifdef _GNU_SOURCE
54#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r)
55#else
56#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp)
57#endif
58
59
60
61extern void pISUpdateComponents(ideal F, const intvec *const V, const int MIN, const ring r);
62// extern ring rCurrRingAssure_SyzComp();
63extern ring rAssure_InducedSchreyerOrdering(const ring r, BOOLEAN complete, int sign);
64extern int rGetISPos(const int p, const ring r);
65
66USING_NAMESPACE( SINGULARXXNAME :: DEBUG )
67USING_NAMESPACE( SINGULARXXNAME :: NF )
68
69
70BEGIN_NAMESPACE_NONAME
71
72
73static inline void NoReturn(leftv& res)
74{
75  res->rtyp = NONE;
76  res->data = NULL;
77}
78
79/// wrapper around n_ClearContent
80static BOOLEAN _ClearContent(leftv res, leftv h)
81{
82  NoReturn(res);
83
84  const char *usage = "'ClearContent' needs a (non-zero!) poly or vector argument...";
85 
86  if( h == NULL )
87  {
88    WarnS(usage);
89    return TRUE;
90  }
91
92  assume( h != NULL );
93
94  if( !( h->Typ() == POLY_CMD || h->Typ() == VECTOR_CMD) )
95  {
96    WarnS(usage);
97    return TRUE;
98  }
99
100  assume (h->Next() == NULL);
101 
102  poly ph = reinterpret_cast<poly>(h->Data());
103 
104  if( ph == NULL )
105  {
106    WarnS(usage);
107    return TRUE;
108  }
109 
110  const ring r =  currRing;
111  assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf;
112
113  number n;
114
115  // experimentall (recursive enumerator treatment) of alg. ext
116  CPolyCoeffsEnumerator itr(ph);
117  n_ClearContent(itr, n, C);
118
119  res->data = n;
120  res->rtyp = NUMBER_CMD;
121
122  return FALSE;
123}
124
125/// wrapper around n_ClearDenominators
126static BOOLEAN _ClearDenominators(leftv res, leftv h)
127{
128  NoReturn(res);
129
130  const char *usage = "'ClearDenominators' needs a (non-zero!) poly or vector argument...";
131
132  if( h == NULL )
133  {
134    WarnS(usage);
135    return TRUE;
136  }
137
138  assume( h != NULL );
139
140  if( !( h->Typ() == POLY_CMD || h->Typ() == VECTOR_CMD) )
141  {
142    WarnS(usage);
143    return TRUE;
144  }
145
146  assume (h->Next() == NULL);
147
148  poly ph = reinterpret_cast<poly>(h->Data());
149
150  if( ph == NULL )
151  {
152    WarnS(usage);
153    return TRUE;
154  }
155
156  const ring r =  currRing;
157  assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf;
158
159  number n;
160
161  // experimentall (recursive enumerator treatment) of alg. ext.
162  CPolyCoeffsEnumerator itr(ph);
163  n_ClearDenominators(itr, n, C);
164
165  res->data = n;
166  res->rtyp = NUMBER_CMD;
167
168  return FALSE;
169}
170
171
172/// try to get an optional (simple) integer argument out of h
173/// or return the default value
174static int getOptionalInteger(const leftv& h, const int _n)
175{
176  if( h!= NULL && h->Typ() == INT_CMD )
177  {
178    int n = (int)(long)(h->Data());
179
180    if( n < 0 )
181      Warn("Negative (%d) optional integer argument", n);
182
183    return (n);
184  }
185
186  return (_n); 
187}
188
189static BOOLEAN noop(leftv __res, leftv /*__v*/)
190{
191  NoReturn(__res);
192  return FALSE;
193}
194
195static BOOLEAN _ProfilerStart(leftv __res, leftv h)
196{
197  NoReturn(__res);
198#if GOOGLE_PROFILE_ENABLED
199  if( h!= NULL && h->Typ() == STRING_CMD )
200  {
201    const char* name = (char*)(h->Data());
202    assume( name != NULL );   
203    ProfilerStart(name);
204  } else
205    WerrorS("ProfilerStart requires a string [name] argument"); 
206#else
207  WarnS("Sorry no google profiler support (GOOGLE_PROFILE_ENABLE!=1)...");
208//  return TRUE; // ?
209#endif // #if GOOGLE_PROFILE_ENABLED
210  return FALSE;
211  (void)h;
212}
213static BOOLEAN _ProfilerStop(leftv __res, leftv /*__v*/)
214{
215  NoReturn(__res);
216#if GOOGLE_PROFILE_ENABLED
217  ProfilerStop();
218#else
219  WarnS("Sorry no google profiler support (GOOGLE_PROFILE_ENABLED!=1)...");
220//  return TRUE; // ?
221#endif // #if GOOGLE_PROFILE_ENABLED
222  return FALSE;
223}
224
225static inline number jjLONG2N(long d)
226{
227  return n_Init(d, coeffs_BIGINT);
228}
229
230static inline void view(const intvec* v)
231{
232#ifndef SING_NDEBUG
233  v->view();
234#else
235  // This code duplication is only due to Hannes's #ifndef SING_NDEBUG!
236  Print ("intvec: {rows: %d, cols: %d, length: %d, Values: \n", v->rows(), v->cols(), v->length());
237
238  for (int i = 0; i < v->rows(); i++)
239  {
240    Print ("Row[%3d]:", i);
241    for (int j = 0; j < v->cols(); j++)
242      Print (" %5d", (*v)[j + i * (v->cols())] );
243    PrintLn ();
244  }
245  PrintS ("}\n");
246#endif
247
248}
249
250                   
251
252static BOOLEAN DetailedPrint(leftv __res, leftv h)
253{
254  NoReturn(__res);
255
256  if( h == NULL )
257  {
258    WarnS("DetailedPrint needs an argument...");
259    return TRUE;
260  }
261
262  if( h->Typ() == NUMBER_CMD)
263  {
264    number n = (number)h->Data(); 
265
266    const ring r = currRing;
267
268#ifdef LDEBUG
269    r->cf->cfDBTest(n,__FILE__,__LINE__,r->cf);
270#endif
271
272    StringSetS("");
273    n_Write(n, r->cf);
274    PrintS(StringEndS());
275    PrintLn();
276
277    return FALSE;
278  }
279 
280  if( h->Typ() == RING_CMD)
281  {
282    const ring r = (const ring)h->Data();
283    rWrite(r, TRUE);
284    PrintLn();
285#ifdef RDEBUG
286    rDebugPrint(r);
287#endif
288    return FALSE;
289  }
290
291  if( h->Typ() == POLY_CMD || h->Typ() == VECTOR_CMD)
292  {
293    const poly p = (const poly)h->Data(); h = h->Next();
294
295    dPrint(p, currRing, currRing, getOptionalInteger(h, 3));
296
297    return FALSE;
298  }
299
300  if( h->Typ() == IDEAL_CMD || h->Typ() == MODUL_CMD)
301  {
302    const ideal id = (const ideal)h->Data(); h = h->Next(); 
303
304    dPrint(id, currRing, currRing, getOptionalInteger(h, 3));
305   
306    return FALSE;           
307  }
308
309  if( h->Typ() == RESOLUTION_CMD )
310  {
311    const syStrategy syzstr = reinterpret_cast<const syStrategy>(h->Data());
312
313    h = h->Next();
314
315    int nTerms = getOptionalInteger(h, 1);
316
317
318    Print("RESOLUTION_CMD(%p): ", reinterpret_cast<const void*>(syzstr)); PrintLn();
319
320    const ring save = currRing;
321    const ring r = syzstr->syRing;
322    const ring rr = (r != NULL) ? r: save;
323
324
325    const int iLength = syzstr->length;
326
327    Print("int 'length': %d", iLength); PrintLn();
328    Print("int 'regularity': %d", syzstr->regularity); PrintLn();
329    Print("short 'list_length': %hd", syzstr->list_length); PrintLn();
330    Print("short 'references': %hd", syzstr->references); PrintLn();
331
332
333#define PRINT_pINTVECTOR(s, v) Print("intvec '%10s'(%p)", #v, reinterpret_cast<const void*>((s)->v)); \
334if( (s)->v != NULL ){ PrintS(": "); view((s)->v); }; \
335PrintLn();
336
337    PRINT_pINTVECTOR(syzstr, resolution);
338    PRINT_pINTVECTOR(syzstr, betti);
339    PRINT_pINTVECTOR(syzstr, Tl);
340    PRINT_pINTVECTOR(syzstr, cw);
341#undef PRINT_pINTVECTOR
342
343    if (r == NULL)
344      Print("ring '%10s': NULL", "syRing");
345    else 
346      if (r == currRing)
347        Print("ring '%10s': currRing", "syRing");
348      else
349        if (r != NULL && r != save)
350        {
351          Print("ring '%10s': ", "syRing");
352          rWrite(r);
353#ifdef RDEBUG
354          //              rDebugPrint(r);
355#endif
356          // rChangeCurrRing(r);
357        }           
358    PrintLn();
359
360    const SRes rP = syzstr->resPairs;
361    Print("SRes 'resPairs': %p", reinterpret_cast<const void*>(rP)); PrintLn();
362
363    if (rP != NULL)
364      for (int iLevel = 0; (iLevel < iLength) && (rP[iLevel] != NULL) && ((*syzstr->Tl)[iLevel] >= 0); iLevel++)
365      {
366        int n = 0;
367        const int iTl = (*syzstr->Tl)[iLevel];
368        for (int j = 0; (j < iTl) && ((rP[iLevel][j].lcm!=NULL) || (rP[iLevel][j].syz!=NULL)); j++)
369        {
370          if (rP[iLevel][j].isNotMinimal==NULL)
371            n++;
372        }
373        Print("minimal-resPairs-Size[1+%d]: %d", iLevel, n); PrintLn();
374      }
375
376
377    //  const ring rrr = (iLevel > 0) ? rr : save; ?
378#define PRINT_RESOLUTION(s, v) Print("resolution '%12s': %p", #v, reinterpret_cast<const void*>((s)->v)); PrintLn(); \
379if ((s)->v != NULL) \
380  for (int iLevel = 0; (iLevel < iLength) && ( ((s)->v)[iLevel] != NULL ); iLevel++) \
381  { \
382    /* const ring rrr = (iLevel > 0) ? save : save; */ \
383    Print("id '%10s'[%d]: (%p) ncols = %d / size: %d; nrows = %d, rank = %ld / rk: %ld", #v, iLevel, reinterpret_cast<const void*>(((s)->v)[iLevel]), ((s)->v)[iLevel]->ncols, idSize(((s)->v)[iLevel]), ((s)->v)[iLevel]->nrows, ((s)->v)[iLevel]->rank, -1L/*id_RankFreeModule(((s)->v)[iLevel], rrr)*/ ); \
384    PrintLn(); \
385  } \
386  PrintLn();
387
388    // resolvente:
389    PRINT_RESOLUTION(syzstr, minres);
390    PRINT_RESOLUTION(syzstr, fullres);
391
392    assume (id_RankFreeModule (syzstr->res[1], rr) == syzstr->res[1]->rank);
393
394    PRINT_RESOLUTION(syzstr, res);
395    PRINT_RESOLUTION(syzstr, orderedRes);
396#undef PRINT_RESOLUTION
397
398#define PRINT_POINTER(s, v) Print("pointer '%17s': %p", #v, reinterpret_cast<const void*>((s)->v)); PrintLn();
399    // 2d arrays:
400    PRINT_POINTER(syzstr, truecomponents);
401    PRINT_POINTER(syzstr, ShiftedComponents);
402    PRINT_POINTER(syzstr, backcomponents);
403    PRINT_POINTER(syzstr, Howmuch);
404    PRINT_POINTER(syzstr, Firstelem);
405    PRINT_POINTER(syzstr, elemLength);
406    PRINT_POINTER(syzstr, sev);
407
408    // arrays of intvects:
409    PRINT_POINTER(syzstr, weights);
410    PRINT_POINTER(syzstr, hilb_coeffs);
411#undef PRINT_POINTER
412
413
414    if (syzstr->fullres==NULL)
415    {
416      PrintS("resolution 'fullres': (NULL) => resolution not computed yet");
417      PrintLn();
418    } else
419    {
420      Print("resolution 'fullres': (%p) => resolution seems to be computed already", reinterpret_cast<const void*>(syzstr->fullres));
421      PrintLn();
422      dPrint(*syzstr->fullres, save, save, nTerms);
423    }
424
425
426
427
428    if (syzstr->minres==NULL)
429    {
430      PrintS("resolution 'minres': (NULL) => resolution not minimized yet");
431      PrintLn();
432    } else
433    {
434      Print("resolution 'minres': (%p) => resolution seems to be minimized already", reinterpret_cast<const void*>(syzstr->minres));
435      PrintLn();
436      dPrint(*syzstr->minres, save, save, nTerms);
437    }
438
439
440
441
442    /*
443    int ** truecomponents;
444    long** ShiftedComponents;
445    int ** backcomponents;
446    int ** Howmuch;
447    int ** Firstelem;
448    int ** elemLength;
449    unsigned long ** sev;
450
451    intvec ** weights;
452    intvec ** hilb_coeffs;
453
454    SRes resPairs;               //polynomial data for internal use only
455
456    resolvente fullres;
457    resolvente minres;
458    resolvente res;              //polynomial data for internal use only
459    resolvente orderedRes;       //polynomial data for internal use only
460*/
461
462    //            if( currRing != save ) rChangeCurrRing(save);
463  }
464
465
466  return FALSE;
467}
468
469
470
471
472
473
474static BOOLEAN Tail(leftv res, leftv h)
475{
476  NoReturn(res);
477
478  if( h == NULL )
479  {
480    WarnS("Tail needs an argument...");
481    return TRUE;
482  }
483
484  assume( h != NULL );
485
486  if( h->Typ() == POLY_CMD || h->Typ() == VECTOR_CMD)
487  {
488    const poly p = (const poly)h->Data();
489
490    res->rtyp = h->Typ();
491
492    if( p == NULL)
493      res->data = NULL;
494    else
495      res->data = p_Copy( pNext(p), currRing );
496
497    h = h->Next(); assume (h == NULL);
498   
499    return FALSE;
500  }
501
502  if( h->Typ() == IDEAL_CMD || h->Typ() == MODUL_CMD)
503  {
504    const ideal id = (const ideal)h->Data();
505
506    res->rtyp = h->Typ();
507
508    if( id == NULL)
509      res->data = NULL;
510    else
511    {
512      const ideal newid = idInit(IDELEMS(id),id->rank);
513      for (int i=IDELEMS(id) - 1; i >= 0; i--)
514        if( id->m[i] != NULL )
515          newid->m[i] = p_Copy(pNext(id->m[i]), currRing);
516        else
517          newid->m[i] = NULL;
518     
519      res->data = newid; 
520    }
521   
522    h = h->Next(); assume (h == NULL);
523   
524    return FALSE;
525  }
526
527  WarnS("Tail needs a single poly/vector/ideal/module argument...");
528  return TRUE;
529}
530
531
532#ifdef _GNU_SOURCE
533static int cmp_c_ds(const void *p1, const void *p2, void *R)
534{
535#else
536static int cmp_c_ds(const void *p1, const void *p2)
537{
538  void *R = currRing; 
539#endif
540
541  const int YES = 1;
542  const int NO = -1;
543
544  const ring r =  (const ring) R; // TODO/NOTE: the structure is known: C, lp!!!
545
546  assume( r == currRing );
547
548  const poly a = *(const poly*)p1;
549  const poly b = *(const poly*)p2;
550
551  assume( a != NULL );
552  assume( b != NULL );
553 
554  assume( p_LmTest(a, r) );
555  assume( p_LmTest(b, r) );
556
557
558  const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r);
559
560  // TODO: test this!!!!!!!!!!!!!!!!
561
562  //return -( compare (c, qsorts) )
563
564  const int __DEBUG__ = 0;
565
566#ifndef NDEBUG
567  if( __DEBUG__ )
568  {
569    PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2);
570    PrintS("b: "); dPrint(b, r, r, 2);
571    PrintLn();
572  }
573#endif
574
575
576  if( iCompDiff > 0 )
577    return YES;
578
579  if( iCompDiff < 0 )
580    return  NO;
581
582  assume( iCompDiff == 0 );
583
584  const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r);
585
586  if( iDegDiff > 0 )
587    return YES;
588
589  if( iDegDiff < 0 )
590    return  NO;
591
592  assume( iDegDiff == 0 );
593
594#ifndef NDEBUG
595  if( __DEBUG__ )
596  {
597    PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn();
598  }
599#endif
600
601  for (int v = rVar(r); v > 0; v--)
602  {
603    assume( v > 0 );
604    assume( v <= rVar(r) );
605
606    const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r);
607
608    if( d > 0 )
609      return YES;
610
611    if( d < 0 )
612      return NO;
613
614    assume( d == 0 );
615  }
616
617  return 0; 
618}
619
620
621
622static ideal ComputeLeadingSyzygyTerms(const ideal& id, const ring r)
623{
624  // 1. set of components S?
625  // 2. for each component c from S: set of indices of leading terms
626  // with this component?
627  // 3. short exp. vectors for each leading term?
628
629  const int size = IDELEMS(id);
630
631  if( size < 2 )
632  {
633    const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal...       
634    return newid;
635  }
636
637
638  // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!??
639
640  // components should come in groups: count elements in each group
641  // && estimate the real size!!!
642
643
644  // use just a vector instead???
645  const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
646
647  int k = 0;
648
649  for (int j = 0; j < size; j++)
650  {
651    const poly p = id->m[j];
652    assume( p != NULL );
653    const int  c = p_GetComp(p, r);
654
655    for (int i = j - 1; i >= 0; i--)
656    {
657      const poly pp = id->m[i];
658      assume( pp != NULL );
659      const int  cc = p_GetComp(pp, r);
660
661      if( c != cc )
662        continue;
663
664      const poly m = p_Init(r); // p_New???
665
666      // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)!
667      for (int v = rVar(r); v > 0; v--)
668      {
669        assume( v > 0 );
670        assume( v <= rVar(r) );
671
672        const short e1 = p_GetExp(p , v, r);
673        const short e2 = p_GetExp(pp, v, r);
674
675        if( e1 >= e2 )
676          p_SetExp(m, v, 0, r);
677        else
678          p_SetExp(m, v, e2 - e1, r);
679
680      }
681
682      assume( (j > i) && (i >= 0) );
683
684      p_SetComp(m, j + 1, r);
685      pNext(m) = NULL;
686      p_SetCoeff0(m, n_Init(1, r->cf), r); // for later...
687
688      p_Setm(m, r); // should not do anything!!!
689
690      newid->m[k++] = m;
691    }
692  }
693
694//   if( __DEBUG__ && FALSE )
695//   {
696//     PrintS("ComputeLeadingSyzygyTerms::Temp0: \n");
697//     dPrint(newid, r, r, 1);
698//   }
699
700  // the rest of newid is assumed to be zeroes...
701
702  // simplify(newid, 2 + 32)??
703  // sort(newid, "C,ds")[1]???
704  id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
705
706//   if( __DEBUG__ && FALSE )
707//   {
708//     PrintS("ComputeLeadingSyzygyTerms::Temp1: \n");
709//     dPrint(newid, r, r, 1);
710//   }
711
712  idSkipZeroes(newid); // #define SIMPL_NULL 2
713
714//   if( __DEBUG__ )
715//   {
716//     PrintS("ComputeLeadingSyzygyTerms::Output: \n");
717//     dPrint(newid, r, r, 1);
718//   }
719
720  const int sizeNew = IDELEMS(newid);
721
722  if( sizeNew >= 2 )
723    qsort_my(newid->m, sizeNew, sizeof(poly), r, cmp_c_ds);
724
725  newid->rank = id_RankFreeModule(newid, r);
726
727  return newid;
728}
729
730
731static BOOLEAN ComputeLeadingSyzygyTerms(leftv res, leftv h)
732{
733  const ring r = currRing;
734  NoReturn(res);
735
736  if( h == NULL )
737  {
738    WarnS("ComputeLeadingSyzygyTerms needs an argument...");
739    return TRUE;
740  }
741
742  assume( h != NULL ); 
743
744#ifndef NDEBUG
745  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
746#else
747  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
748#endif
749
750  if( h->Typ() == IDEAL_CMD || h->Typ() == MODUL_CMD)
751  {
752    const ideal id = (const ideal)h->Data();
753
754    assume(id != NULL);
755
756    if( __DEBUG__ )
757    {
758      PrintS("ComputeLeadingSyzygyTerms::Input: \n");
759     
760      const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)0)));
761      const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
762      const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));
763
764
765      Print("\nSYZCHECK: \t%d", __SYZCHECK__);
766      Print(", DEBUG: \t%d", __DEBUG__);
767      Print(", LEAD2SYZ: \t%d", __LEAD2SYZ__);
768      Print(", TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
769
770      dPrint(id, r, r, 1);
771
772      assume( !__LEAD2SYZ__ );
773    }
774
775    h = h->Next(); assume (h == NULL);
776
777    const ideal newid = ComputeLeadingSyzygyTerms(id, r);
778   
779    res->data = newid; res->rtyp = MODUL_CMD;
780    return FALSE;
781  }
782
783  WarnS("ComputeLeadingSyzygyTerms needs a single ideal/module argument...");
784  return TRUE;
785}
786
787///  sorting wrt <c,ds> & reversing...
788/// change the input inplace!!!
789// TODO: use a ring with >_{c, ds}!???
790static BOOLEAN Sort_c_ds(leftv res, leftv h)
791{
792  NoReturn(res);
793
794  const ring r = currRing;
795  NoReturn(res);
796
797  if( h == NULL )
798  {
799    WarnS("Sort_c_ds needs an argument...");
800    return TRUE;
801  }
802
803  assume( h != NULL ); 
804
805#ifndef NDEBUG
806  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
807#else
808  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
809#endif
810
811  if(    (h->Typ() == IDEAL_CMD || h->Typ() == MODUL_CMD)
812      && (h->rtyp  == IDHDL) // must be a variable!
813      && (h->e == NULL) // not a list element
814      ) 
815  {
816    const ideal id = (const ideal)h->Data();
817
818    assume(id != NULL);
819
820    if( __DEBUG__ )
821    {
822      PrintS("Sort_c_ds::Input: \n");
823
824      const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)0)));
825      const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
826      const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));   
827
828      Print("\nSYZCHECK: \t%d", __SYZCHECK__);
829      Print(", DEBUG: \t%d", __DEBUG__);
830      Print(", LEAD2SYZ: \t%d", __LEAD2SYZ__);
831      Print(", TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
832
833      dPrint(id, r, r, 1);     
834    }
835
836    assume (h->Next() == NULL);
837
838    id_Test(id, r);
839
840    const int size = IDELEMS(id);
841
842    const ideal newid = id; // id_Copy(id, r); // copy???
843
844    if( size >= 2 )
845      qsort_my(newid->m, size, sizeof(poly), r, cmp_c_ds);
846   
847//    res->data = newid;
848//    res->rtyp = h->Typ();
849   
850    if( __DEBUG__ )
851    {
852      PrintS("Sort_c_ds::Output: \n");
853      dPrint(newid, r, r, 1);
854    }
855
856    return FALSE;
857  }
858
859  WarnS("ComputeLeadingSyzygyTerms needs a single ideal/module argument (must be a variable!)...");
860  return TRUE; 
861}
862
863
864
865static ideal Compute2LeadingSyzygyTerms(const ideal& id, const ring r)
866{
867#ifndef NDEBUG
868  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
869#else
870  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
871#endif
872
873  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
874  const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));   
875
876
877
878  // 1. set of components S?
879  // 2. for each component c from S: set of indices of leading terms
880  // with this component?
881  // 3. short exp. vectors for each leading term?
882
883  const int size = IDELEMS(id);
884
885  if( size < 2 )
886  {
887    const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...   
888    return newid;
889  }
890
891
892    // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
893
894    // components should come in groups: count elements in each group
895    // && estimate the real size!!!
896
897
898    // use just a vector instead???
899  ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
900
901  int k = 0;
902
903  for (int j = 0; j < size; j++)
904  {
905    const poly p = id->m[j];
906    assume( p != NULL );
907    const int  c = p_GetComp(p, r);
908
909    for (int i = j - 1; i >= 0; i--)
910    {
911      const poly pp = id->m[i];
912      assume( pp != NULL );
913      const int  cc = p_GetComp(pp, r);
914
915      if( c != cc )
916        continue;
917
918        // allocate memory & zero it out!
919      const poly m = p_Init(r); const poly mm = p_Init(r);
920
921
922        // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
923        // TODO: optimize: knowing the ring structure: (C/lp)!
924
925      for (int v = rVar(r); v > 0; v--)
926      {
927        assume( v > 0 );
928        assume( v <= rVar(r) );
929
930        const short e1 = p_GetExp(p , v, r);
931        const short e2 = p_GetExp(pp, v, r);
932
933        if( e1 >= e2 )
934          p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
935        else
936          p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
937
938      }
939
940      assume( (j > i) && (i >= 0) );
941
942      p_SetComp(m, j + 1, r);
943      p_SetComp(mm, i + 1, r);
944
945      const number& lc1 = p_GetCoeff(p , r);
946      const number& lc2 = p_GetCoeff(pp, r);
947
948      number g = n_Lcm( lc1, lc2, r );
949
950      p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
951      p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);
952
953      n_Delete(&g, r);
954
955      p_Setm(m, r); // should not do anything!!!
956      p_Setm(mm, r); // should not do anything!!!
957
958      pNext(m) = mm; //        pNext(mm) = NULL;
959
960      newid->m[k++] = m;
961    }
962  }
963
964//   if( __DEBUG__ && FALSE )
965//   {
966//     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
967//     dPrint(newid, r, r, 1);
968//   }
969
970  if( !__TAILREDSYZ__ )
971  {
972      // simplify(newid, 2 + 32)??
973      // sort(newid, "C,ds")[1]???
974    id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
975
976//     if( __DEBUG__ && FALSE )
977//     {
978//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
979//       dPrint(newid, r, r, 1);
980//     }
981  }
982  else
983  {
984      //      option(redSB); option(redTail);
985      //      TEST_OPT_REDSB
986      //      TEST_OPT_REDTAIL
987    assume( r == currRing );
988    BITSET _save_test = test;
989    test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB));
990
991    intvec* w=new intvec(IDELEMS(newid));
992    ideal tmp = kStd(newid, currQuotient, isHomog, &w);
993    delete w;
994
995    test = _save_test;
996
997    id_Delete(&newid, r);
998    newid = tmp;
999
1000//     if( __DEBUG__ && FALSE )
1001//     {
1002//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
1003//       dPrint(newid, r, r, 1);
1004//     }
1005
1006  }
1007
1008  idSkipZeroes(newid);
1009
1010  const int sizeNew = IDELEMS(newid);
1011
1012  if( sizeNew >= 2 )
1013    qsort_my(newid->m, sizeNew, sizeof(poly), r, cmp_c_ds);
1014
1015  newid->rank = id_RankFreeModule(newid, r);
1016 
1017  return newid;
1018}
1019
1020
1021
1022static BOOLEAN Compute2LeadingSyzygyTerms(leftv res, leftv h)
1023{
1024  const ring r = currRing;
1025  NoReturn(res);
1026
1027  if( h == NULL )
1028  {
1029    WarnS("Compute2LeadingSyzygyTerms needs an argument...");
1030    return TRUE;
1031  }
1032
1033  assume( h != NULL ); 
1034
1035#ifndef NDEBUG
1036  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
1037#else
1038  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
1039#endif
1040
1041  const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)0)));
1042 
1043  assume( __LEAD2SYZ__ );
1044 
1045  if( h->Typ() == IDEAL_CMD || h->Typ() == MODUL_CMD)
1046  {
1047    const ideal id = (const ideal)h->Data();
1048
1049    assume(id != NULL);
1050
1051    if( __DEBUG__ )
1052    {
1053      PrintS("Compute2LeadingSyzygyTerms::Input: \n");
1054      dPrint(id, r, r, 0);
1055    }
1056
1057    h = h->Next(); assume (h == NULL);
1058
1059    res->data = Compute2LeadingSyzygyTerms(id, r);
1060    res->rtyp = MODUL_CMD;
1061
1062    return FALSE;
1063  }
1064
1065  WarnS("Compute2LeadingSyzygyTerms needs a single ideal/module argument...");
1066  return TRUE;
1067}
1068
1069
1070/// return a new term: leading coeff * leading monomial of p
1071/// with 0 leading component!
1072static inline poly leadmonom(const poly p, const ring r)
1073{
1074  poly m = NULL;
1075
1076  if( p != NULL )
1077  {
1078    assume( p != NULL );
1079    assume( p_LmTest(p, r) );
1080
1081    m = p_LmInit(p, r);
1082    p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);
1083
1084    p_SetComp(m, 0, r);
1085    p_Setm(m, r);
1086
1087    assume( p_GetComp(m, r) == 0 );
1088    assume( m != NULL );
1089    assume( pNext(m) == NULL );
1090    assume( p_LmTest(m, r) );
1091  }
1092
1093  return m;
1094}
1095
1096
1097
1098static poly FindReducer(poly product, poly syzterm,
1099                        ideal L, ideal LS,
1100                        const ring r)
1101{
1102  assume( product != NULL );
1103  assume( L != NULL );
1104
1105  int c = 0;
1106
1107  if (syzterm != NULL)
1108    c = p_GetComp(syzterm, r) - 1;
1109
1110  assume( c >= 0 && c < IDELEMS(L) );
1111 
1112#ifndef NDEBUG
1113  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
1114#else
1115  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
1116#endif
1117
1118  long debug = __DEBUG__;
1119  const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)debug)));   
1120 
1121  if (__SYZCHECK__ && syzterm != NULL)
1122  {
1123    const poly m = L->m[c];
1124
1125    assume( m != NULL ); assume( pNext(m) == NULL );
1126
1127    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
1128    assume( p_EqualPolys(lm, product, r) );
1129
1130    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
1131    //  def @@product = leadmonomial(syzterm) * L[@@r];
1132
1133    p_Delete(&lm, r);   
1134  }
1135 
1136  // looking for an appropriate diviser q = L[k]...
1137  for( int k = IDELEMS(L)-1; k>= 0; k-- )
1138  {
1139    const poly p = L->m[k];   
1140
1141    // ... which divides the product, looking for the _1st_ appropriate one!
1142    if( !p_LmDivisibleBy(p, product, r) )
1143      continue;
1144
1145
1146    const poly q = p_New(r);
1147    pNext(q) = NULL;
1148    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
1149    p_SetComp(q, k + 1, r);
1150    p_Setm(q, r);
1151
1152    // cannot allow something like: a*gen(i) - a*gen(i)
1153    if (syzterm != NULL && (k == c))
1154      if (p_ExpVectorEqual(syzterm, q, r))
1155      {
1156        if( __DEBUG__ )
1157        {
1158          Print("FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1159          dPrint(syzterm, r, r, 1);
1160        }   
1161
1162        p_LmFree(q, r);
1163        continue;
1164      }
1165
1166    // while the complement (the fraction) is not reducible by leading syzygies
1167    if( LS != NULL )
1168    {
1169      BOOLEAN ok = TRUE;
1170
1171      for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )
1172      {
1173        const poly pp = LS->m[kk];
1174
1175        if( p_LmDivisibleBy(pp, q, r) )
1176        {
1177         
1178          if( __DEBUG__ )
1179          {
1180            Print("FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);
1181            dPrint(pp, r, r, 1);
1182          }   
1183
1184          ok = FALSE; // q in <LS> :((
1185          break;                 
1186        }
1187      }
1188
1189      if(!ok)
1190      {
1191        p_LmFree(q, r);
1192        continue;
1193      }
1194    }
1195
1196    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
1197    return q;
1198
1199  }
1200
1201 
1202  return NULL;
1203}
1204
1205
1206/// TODO: save shortcut (syz: |-.->) LM(LM(m) * "t") -> syz?
1207/// proc SSFindReducer(def product, def syzterm, def L, def T, list #)
1208static BOOLEAN FindReducer(leftv res, leftv h)
1209{
1210  const char* usage = "`FindReducer(<poly/vector>, <vector/0>, <ideal/module>[,<module>])` expected";
1211  const ring r = currRing;
1212
1213  NoReturn(res);
1214
1215
1216#ifndef NDEBUG
1217  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
1218#else
1219  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
1220#endif
1221
1222  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
1223
1224  if ((h==NULL) || (h->Typ()!=VECTOR_CMD && h->Typ() !=POLY_CMD) || (h->Data() == NULL))
1225  {
1226    WerrorS(usage);
1227    return TRUE;   
1228  }
1229   
1230  const poly product = (poly) h->Data(); assume (product != NULL);
1231
1232
1233  h = h->Next();
1234  if ((h==NULL) || !((h->Typ()==VECTOR_CMD) || (h->Data() == NULL)) )
1235  {
1236    WerrorS(usage);
1237    return TRUE;   
1238  }
1239
1240  poly syzterm = NULL;
1241
1242  if(h->Typ()==VECTOR_CMD) 
1243    syzterm = (poly) h->Data();
1244
1245
1246
1247  h = h->Next();
1248  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1249  {
1250    WerrorS(usage);
1251    return TRUE;   
1252  }
1253 
1254  const ideal L = (ideal) h->Data(); h = h->Next();
1255
1256  assume( IDELEMS(L) > 0 );
1257
1258  ideal LS = NULL;
1259
1260  if ((h != NULL) && (h->Typ() ==MODUL_CMD) && (h->Data() != NULL))
1261  {
1262    LS = (ideal)h->Data();
1263    h = h->Next();
1264  }
1265
1266  if( __TAILREDSYZ__ )
1267    assume (LS != NULL);
1268
1269  assume( h == NULL );
1270
1271  if( __DEBUG__ )
1272  {
1273    PrintS("FindReducer(product, syzterm, L, T, #)::Input: \n");
1274
1275    PrintS("product: "); dPrint(product, r, r, 2);
1276    PrintS("syzterm: "); dPrint(syzterm, r, r, 2);
1277    PrintS("L: "); dPrint(L, r, r, 0);
1278//    PrintS("T: "); dPrint(T, r, r, 4);
1279
1280    if( LS == NULL )
1281      PrintS("LS: NULL\n");
1282    else
1283    {
1284      PrintS("LS: "); dPrint(LS, r, r, 0);
1285    }
1286  }
1287
1288  res->rtyp = VECTOR_CMD;
1289  res->data = FindReducer(product, syzterm, L, LS, r);
1290
1291  if( __DEBUG__ )
1292  {
1293    PrintS("FindReducer::Output: \n");
1294    dPrint((poly)res->data, r, r, 2);
1295  }   
1296 
1297  return FALSE;   
1298 
1299}
1300
1301static poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
1302                       ideal L, ideal T, ideal LS, const ring r);
1303
1304static poly TraverseTail(poly multiplier, poly tail, 
1305                         ideal L, ideal T, ideal LS,
1306                         const ring r)
1307{
1308  assume( multiplier != NULL );
1309
1310  assume( L != NULL );
1311  assume( T != NULL );
1312
1313  poly s = NULL;
1314
1315  // iterate over the tail
1316  for(poly p = tail; p != NULL; p = pNext(p))
1317    s = p_Add_q(s, ReduceTerm(multiplier, p, NULL, L, T, LS, r), r); 
1318   
1319  return s;
1320}
1321
1322
1323static poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
1324                       ideal L, ideal T, ideal LS, const ring r)
1325{
1326  assume( multiplier != NULL );
1327  assume( term4reduction != NULL );
1328
1329
1330  assume( L != NULL );
1331  assume( T != NULL );
1332 
1333  // assume(r == currRing); // ?
1334
1335  // simple implementation with FindReducer:
1336  poly s = NULL;
1337 
1338  if(1)
1339  {
1340    // NOTE: only LT(term4reduction) should be used in the following:
1341    poly product = pp_Mult_mm(multiplier, term4reduction, r);
1342    s = FindReducer(product, syztermCheck, L, LS, r);
1343    p_Delete(&product, r);
1344  }
1345
1346  if( s == NULL ) // No Reducer?
1347    return s;
1348
1349  poly b = leadmonom(s, r);
1350
1351  const int c = p_GetComp(s, r) - 1;
1352  assume( c >= 0 && c < IDELEMS(T) );
1353
1354  const poly tail = T->m[c];
1355
1356  if( tail != NULL )
1357    s = p_Add_q(s, TraverseTail(b, tail, L, T, LS, r), r); 
1358 
1359  return s;
1360}
1361
1362
1363static poly SchreyerSyzygyNF(poly syz_lead, poly syz_2, ideal L, ideal T, ideal LS, const ring r)
1364{
1365  assume( syz_lead != NULL );
1366  assume( syz_2 != NULL );
1367
1368  assume( L != NULL );
1369  assume( T != NULL );
1370
1371  assume( IDELEMS(L) == IDELEMS(T) );
1372
1373  int  c = p_GetComp(syz_lead, r) - 1;
1374
1375  assume( c >= 0 && c < IDELEMS(T) );
1376
1377  poly p = leadmonom(syz_lead, r); // :( 
1378  poly spoly = pp_Mult_qq(p, T->m[c], r);
1379  p_Delete(&p, r);
1380
1381 
1382  c = p_GetComp(syz_2, r) - 1;
1383  assume( c >= 0 && c < IDELEMS(T) );
1384
1385  p = leadmonom(syz_2, r); // :(
1386  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
1387  p_Delete(&p, r);
1388
1389  poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
1390
1391  while (spoly != NULL)
1392  {
1393    poly t = FindReducer(spoly, NULL, L, LS, r);
1394
1395    p_LmDelete(&spoly, r);
1396   
1397    if( t != NULL )
1398    {
1399      p = leadmonom(t, r); // :(
1400      c = p_GetComp(t, r) - 1;
1401
1402      assume( c >= 0 && c < IDELEMS(T) );
1403     
1404      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
1405
1406      p_Delete(&p, r);
1407     
1408      tail = p_Add_q(tail, t, r);
1409    }   
1410  }
1411 
1412  return tail;
1413}
1414
1415// proc SchreyerSyzygyNF(vector syz_lead, vector syz_2, def L, def T, list #)
1416static BOOLEAN SchreyerSyzygyNF(leftv res, leftv h)
1417{
1418  const char* usage = "`SchreyerSyzygyNF(<vector>, <vector>, <ideal/module>, <ideal/module>[,<module>])` expected";
1419  const ring r = currRing;
1420
1421  NoReturn(res);
1422
1423#ifndef NDEBUG
1424  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
1425#else
1426  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
1427#endif
1428
1429  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
1430//  const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)0)));
1431  const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));   
1432
1433  if ((h==NULL) || (h->Typ() != VECTOR_CMD) || (h->Data() == NULL))
1434  {
1435    WerrorS(usage);
1436    return TRUE;   
1437  }
1438
1439  const poly syz_lead = (poly) h->Data(); assume (syz_lead != NULL);
1440
1441
1442  h = h->Next();
1443  if ((h==NULL) || (h->Typ() != VECTOR_CMD) || (h->Data() == NULL))
1444  {
1445    WerrorS(usage);
1446    return TRUE;   
1447  }
1448
1449  const poly syz_2 = (poly) h->Data(); assume (syz_2 != NULL);
1450
1451  h = h->Next();
1452  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1453  {
1454    WerrorS(usage);
1455    return TRUE;   
1456  }
1457
1458  const ideal L = (ideal) h->Data(); assume( IDELEMS(L) > 0 );
1459
1460
1461  h = h->Next();
1462  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1463  {
1464    WerrorS(usage);
1465    return TRUE;   
1466  }
1467
1468  const ideal T = (ideal) h->Data();
1469
1470  assume( IDELEMS(L) == IDELEMS(T) );
1471
1472  ideal LS = NULL;
1473
1474  h = h->Next();
1475  if ((h != NULL) && (h->Typ() ==MODUL_CMD) && (h->Data() != NULL))
1476  {
1477    LS = (ideal)h->Data();
1478    h = h->Next();
1479  }
1480
1481  if( __TAILREDSYZ__ )
1482    assume (LS != NULL);
1483
1484  assume( h == NULL );
1485
1486  if( __DEBUG__ )
1487  {
1488    PrintS("SchreyerSyzygyNF(syz_lead, syz_2, L, T, #)::Input: \n");
1489
1490    PrintS("syz_lead: "); dPrint(syz_lead, r, r, 2);
1491    PrintS("syz_2: "); dPrint(syz_2, r, r, 2);
1492
1493    PrintS("L: "); dPrint(L, r, r, 0);
1494    PrintS("T: "); dPrint(T, r, r, 0);
1495
1496    if( LS == NULL )
1497      PrintS("LS: NULL\n");
1498    else
1499    {
1500      PrintS("LS: "); dPrint(LS, r, r, 0);
1501    }
1502  }
1503 
1504  res->rtyp = VECTOR_CMD;
1505  res->data = SchreyerSyzygyNF(syz_lead, syz_2, L, T, LS, r);
1506
1507  if( __DEBUG__ )
1508  {
1509    PrintS("SchreyerSyzygyNF::Output: ");
1510
1511    dPrint((poly)res->data, r, r, 2);
1512  }
1513
1514
1515  return FALSE;
1516}
1517
1518
1519
1520/// TODO: save shortcut (syz: |-.->) LM(m) * "t" -> ?
1521/// proc SSReduceTerm(poly m, def t, def syzterm, def L, def T, list #)
1522static BOOLEAN ReduceTerm(leftv res, leftv h)
1523{
1524  const char* usage = "`ReduceTerm(<poly>, <poly/vector>, <vector/0>, <ideal/module>, <ideal/module>[,<module>])` expected";
1525  const ring r = currRing;
1526
1527  NoReturn(res);
1528
1529#ifndef NDEBUG
1530  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
1531#else
1532  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
1533#endif
1534
1535  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
1536//  const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)0)));
1537  const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));   
1538
1539  if ((h==NULL) || (h->Typ() !=POLY_CMD) || (h->Data() == NULL))
1540  {
1541    WerrorS(usage);
1542    return TRUE;   
1543  }
1544
1545  const poly multiplier = (poly) h->Data(); assume (multiplier != NULL);
1546
1547 
1548  h = h->Next();
1549  if ((h==NULL) || (h->Typ()!=VECTOR_CMD && h->Typ() !=POLY_CMD) || (h->Data() == NULL))
1550  {
1551    WerrorS(usage);
1552    return TRUE;   
1553  }
1554
1555  const poly term4reduction = (poly) h->Data(); assume( term4reduction != NULL );
1556
1557 
1558  poly syztermCheck = NULL;
1559 
1560  h = h->Next();
1561  if ((h==NULL) || !((h->Typ()==VECTOR_CMD) || (h->Data() == NULL)) )
1562  {
1563    WerrorS(usage);
1564    return TRUE;   
1565  }
1566
1567  if(h->Typ()==VECTOR_CMD) 
1568    syztermCheck = (poly) h->Data();
1569
1570 
1571  h = h->Next();
1572  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1573  {
1574    WerrorS(usage);
1575    return TRUE;   
1576  }
1577
1578  const ideal L = (ideal) h->Data(); assume( IDELEMS(L) > 0 );
1579
1580 
1581  h = h->Next();
1582  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1583  {
1584    WerrorS(usage);
1585    return TRUE;   
1586  }
1587
1588  const ideal T = (ideal) h->Data();
1589
1590  assume( IDELEMS(L) == IDELEMS(T) );
1591
1592  ideal LS = NULL;
1593
1594  h = h->Next();
1595  if ((h != NULL) && (h->Typ() ==MODUL_CMD) && (h->Data() != NULL))
1596  {
1597    LS = (ideal)h->Data();
1598    h = h->Next();
1599  }
1600
1601  if( __TAILREDSYZ__ )
1602    assume (LS != NULL);
1603
1604  assume( h == NULL );
1605
1606  if( __DEBUG__ )
1607  {
1608    PrintS("ReduceTerm(m, t, syzterm, L, T, #)::Input: \n");
1609
1610    PrintS("m: "); dPrint(multiplier, r, r, 2);
1611    PrintS("t: "); dPrint(term4reduction, r, r, 2);
1612    PrintS("syzterm: "); dPrint(syztermCheck, r, r, 2);
1613   
1614    PrintS("L: "); dPrint(L, r, r, 0);
1615    PrintS("T: "); dPrint(T, r, r, 0);
1616
1617    if( LS == NULL )
1618      PrintS("LS: NULL\n");
1619    else
1620    {
1621      PrintS("LS: "); dPrint(LS, r, r, 0);
1622    }
1623  }
1624
1625
1626  if (__SYZCHECK__ && syztermCheck != NULL)
1627  {
1628    const int c = p_GetComp(syztermCheck, r) - 1;
1629    assume( c >= 0 && c < IDELEMS(L) );
1630   
1631    const poly p = L->m[c];
1632    assume( p != NULL ); assume( pNext(p) == NULL );   
1633
1634    assume( p_EqualPolys(term4reduction, p, r) ); // assume? TODO
1635
1636
1637    poly m = leadmonom(syztermCheck, r);
1638    assume( m != NULL ); assume( pNext(m) == NULL );
1639
1640    assume( p_EqualPolys(multiplier, m, r) ); // assume? TODO
1641
1642    p_Delete(&m, r);   
1643   
1644// NOTE:   leadmonomial(syzterm) == m &&  L[leadcomp(syzterm)] == t
1645  }
1646
1647  res->rtyp = VECTOR_CMD;
1648  res->data = ReduceTerm(multiplier, term4reduction, syztermCheck, L, T, LS, r);
1649
1650
1651  if( __DEBUG__ )
1652  {
1653    PrintS("ReduceTerm::Output: ");
1654
1655    dPrint((poly)res->data, r, r, 2);
1656  }
1657 
1658 
1659  return FALSE;
1660}
1661
1662
1663
1664
1665// TODO: store m * @tail -.-^-.-^-.--> ?
1666// proc SSTraverseTail(poly m, def @tail, def L, def T, list #)
1667static BOOLEAN TraverseTail(leftv res, leftv h)
1668{
1669  const char* usage = "`TraverseTail(<poly>, <poly/vector>, <ideal/module>, <ideal/module>[,<module>])` expected";
1670  const ring r = currRing;
1671
1672  NoReturn(res);
1673
1674#ifndef NDEBUG
1675  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
1676#else
1677  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
1678#endif
1679
1680  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)1)));
1681
1682  if ((h==NULL) || (h->Typ() !=POLY_CMD) || (h->Data() == NULL))
1683  {
1684    WerrorS(usage);
1685    return TRUE;   
1686  }
1687
1688  const poly multiplier = (poly) h->Data(); assume (multiplier != NULL);
1689
1690  h = h->Next();
1691  if ((h==NULL) || (h->Typ()!=VECTOR_CMD && h->Typ() !=POLY_CMD))
1692  {
1693    WerrorS(usage);
1694    return TRUE;   
1695  }
1696
1697  const poly tail = (poly) h->Data(); 
1698
1699  h = h->Next();
1700
1701  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1702  {
1703    WerrorS(usage);
1704    return TRUE;   
1705  }
1706
1707  const ideal L = (ideal) h->Data();
1708
1709  assume( IDELEMS(L) > 0 );
1710
1711  h = h->Next();
1712  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1713  {
1714    WerrorS(usage);
1715    return TRUE;   
1716  }
1717
1718  const ideal T = (ideal) h->Data();
1719
1720  assume( IDELEMS(L) == IDELEMS(T) );
1721
1722  h = h->Next();
1723 
1724  ideal LS = NULL;
1725
1726  if ((h != NULL) && (h->Typ() ==MODUL_CMD) && (h->Data() != NULL))
1727  {
1728    LS = (ideal)h->Data();
1729    h = h->Next();
1730  }
1731
1732  if( __TAILREDSYZ__ )
1733    assume (LS != NULL);
1734
1735  assume( h == NULL );
1736
1737  if( __DEBUG__ )
1738  {
1739    PrintS("TraverseTail(m, t, L, T, #)::Input: \n");
1740
1741    PrintS("m: "); dPrint(multiplier, r, r, 2);
1742    PrintS("t: "); dPrint(tail, r, r, 10);
1743
1744    PrintS("L: "); dPrint(L, r, r, 0);
1745    PrintS("T: "); dPrint(T, r, r, 0);
1746
1747    if( LS == NULL )
1748      PrintS("LS: NULL\n");
1749    else
1750    {
1751      PrintS("LS: "); dPrint(LS, r, r, 0);
1752    }
1753  }
1754
1755  res->rtyp = VECTOR_CMD;
1756  res->data = TraverseTail(multiplier, tail, L, T, LS, r);
1757
1758
1759  if( __DEBUG__ )
1760  {
1761    PrintS("TraverseTail::Output: ");
1762    dPrint((poly)res->data, r, r, 2);
1763  }
1764
1765  return FALSE;
1766}
1767
1768
1769
1770
1771
1772
1773/// Get leading term without a module component
1774static BOOLEAN leadmonom(leftv res, leftv h)
1775{
1776  NoReturn(res);
1777
1778  if ((h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD) && (h->Data() != NULL))
1779  {
1780    const ring r = currRing;
1781    const poly p = (poly)(h->Data());
1782
1783    const poly m = leadmonom(p, r);
1784
1785    res->rtyp = POLY_CMD;
1786    res->data = reinterpret_cast<void *>(m);
1787
1788    return FALSE;
1789  }
1790
1791  WerrorS("`leadmonom(<poly/vector>)` expected");
1792  return TRUE;
1793}
1794
1795/// Get leading component
1796static BOOLEAN leadcomp(leftv res, leftv h)
1797{
1798  NoReturn(res);
1799
1800  if ((h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD))
1801  {
1802    const ring r = currRing;
1803
1804    const poly p = (poly)(h->Data());
1805
1806    if (p != NULL )
1807    {
1808      assume( p != NULL );
1809      assume( p_LmTest(p, r) );
1810
1811      const unsigned long iComp = p_GetComp(p, r);
1812
1813  //    assume( iComp > 0 ); // p is a vector
1814
1815      res->data = reinterpret_cast<void *>(jjLONG2N(iComp));
1816    } else
1817      res->data = reinterpret_cast<void *>(jjLONG2N(0));
1818     
1819
1820    res->rtyp = BIGINT_CMD;
1821    return FALSE;
1822  }
1823
1824  WerrorS("`leadcomp(<poly/vector>)` expected");
1825  return TRUE;
1826}
1827
1828
1829
1830
1831/// Get raw leading exponent vector
1832static BOOLEAN leadrawexp(leftv res, leftv h)
1833{
1834  NoReturn(res);
1835
1836  if ((h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD) && (h->Data() != NULL))
1837  {
1838    const ring r = currRing;
1839    const poly p = (poly)(h->Data());
1840
1841    assume( p != NULL );
1842    assume( p_LmTest(p, r) );
1843
1844    const int iExpSize = r->ExpL_Size;
1845
1846//    intvec *iv = new intvec(iExpSize);
1847
1848    lists l=(lists)omAllocBin(slists_bin);
1849    l->Init(iExpSize);
1850
1851    for(int i = iExpSize-1; i >= 0; i--)
1852    {
1853      l->m[i].rtyp = BIGINT_CMD;
1854      l->m[i].data = reinterpret_cast<void *>(jjLONG2N(p->exp[i])); // longs...
1855    }
1856
1857    res->rtyp = LIST_CMD; // list of bigints
1858    res->data = reinterpret_cast<void *>(l);
1859    return FALSE;
1860  }
1861
1862  WerrorS("`leadrawexp(<poly/vector>)` expected");
1863  return TRUE;
1864}
1865
1866
1867/// Endowe the current ring with additional (leading) Syz-component ordering
1868static BOOLEAN MakeSyzCompOrdering(leftv res, leftv /*h*/)
1869{
1870
1871  NoReturn(res);
1872
1873  //    res->data = rCurrRingAssure_SyzComp(); // changes current ring! :(
1874  res->data = reinterpret_cast<void *>(rAssure_SyzComp(currRing, TRUE));
1875  res->rtyp = RING_CMD; // return new ring!
1876  // QRING_CMD?
1877
1878  return FALSE;
1879}
1880
1881
1882/// Same for Induced Schreyer ordering (ordering on components is defined by sign!)
1883static BOOLEAN MakeInducedSchreyerOrdering(leftv res, leftv h)
1884{
1885
1886  NoReturn(res);
1887
1888  int sign = 1;
1889  if ((h!=NULL) && (h->Typ()==INT_CMD))
1890  {
1891    const int s = (int)((long)(h->Data()));
1892
1893    if( s != -1 && s != 1 )
1894    {
1895      WerrorS("`MakeInducedSchreyerOrdering(<int>)` called with wrong integer argument (must be +-1)!");
1896      return TRUE;
1897    }
1898
1899    sign = s;           
1900  }
1901
1902  assume( sign == 1 || sign == -1 );
1903  res->data = reinterpret_cast<void *>(rAssure_InducedSchreyerOrdering(currRing, TRUE, sign));
1904  res->rtyp = RING_CMD; // return new ring!
1905  // QRING_CMD?
1906  return FALSE;
1907}
1908
1909
1910/// Returns old SyzCompLimit, can set new limit
1911static BOOLEAN SetSyzComp(leftv res, leftv h)
1912{
1913  NoReturn(res);
1914
1915  const ring r = currRing;
1916
1917  if( !rIsSyzIndexRing(r) )
1918  {
1919    WerrorS("`SetSyzComp(<int>)` called on incompatible ring (not created by 'MakeSyzCompOrdering'!)");
1920    return TRUE;
1921  }
1922
1923  res->rtyp = INT_CMD;
1924  res->data = reinterpret_cast<void *>(rGetCurrSyzLimit(r)); // return old syz limit
1925
1926  if ((h!=NULL) && (h->Typ()==INT_CMD))
1927  {
1928    const int iSyzComp = (int)reinterpret_cast<long>(h->Data());
1929    assume( iSyzComp > 0 );
1930    rSetSyzComp(iSyzComp, currRing);
1931  }
1932
1933  return FALSE;
1934}
1935
1936/// ?
1937static BOOLEAN GetInducedData(leftv res, leftv h)
1938{
1939  NoReturn(res);
1940
1941  const ring r = currRing;
1942
1943  int p = 0; // which IS-block? p^th!
1944
1945  if ((h!=NULL) && (h->Typ()==INT_CMD))
1946  {
1947    p = (int)((long)(h->Data())); h=h->next;
1948    assume(p >= 0);
1949  }
1950
1951  const int pos = rGetISPos(p, r);
1952
1953  if(  /*(*/ -1 == pos /*)*/  )
1954  {
1955    WerrorS("`GetInducedData([int])` called on incompatible ring (not created by 'MakeInducedSchreyerOrdering'!)");
1956    return TRUE;
1957  }
1958
1959
1960  const int iLimit = r->typ[pos].data.is.limit;
1961  const ideal F = r->typ[pos].data.is.F;
1962  ideal FF = id_Copy(F, r);
1963
1964  lists l=(lists)omAllocBin(slists_bin);
1965  l->Init(2);
1966
1967  l->m[0].rtyp = INT_CMD;
1968  l->m[0].data = reinterpret_cast<void *>(iLimit);
1969
1970
1971  //        l->m[1].rtyp = MODUL_CMD;
1972
1973  if( idIsModule(FF, r) )
1974  {
1975    l->m[1].rtyp = MODUL_CMD;
1976
1977    //          Print("before: %d\n", FF->nrows);
1978    //          FF->nrows = id_RankFreeModule(FF, r); // ???
1979    //          Print("after: %d\n", FF->nrows);
1980  }
1981  else
1982    l->m[1].rtyp = IDEAL_CMD;
1983
1984  l->m[1].data = reinterpret_cast<void *>(FF);
1985
1986  res->rtyp = LIST_CMD; // list of int/module
1987  res->data = reinterpret_cast<void *>(l);
1988
1989  return FALSE;
1990
1991}
1992
1993
1994/* // the following turned out to be unnecessary...   
1995/// Finds p^th AM ordering, and returns its position in r->typ[] AND
1996/// corresponding &r->wvhdl[]
1997/// returns FALSE if something went wrong!
1998/// p - starts with 0!
1999BOOLEAN rGetAMPos(const ring r, const int p, int &typ_pos, int &wvhdl_pos, const BOOLEAN bSearchWvhdl = FALSE)
2000{
2001#if MYTEST
2002  Print("rGetAMPos(p: %d...)\nF:", p);
2003  PrintLn();
2004#endif
2005  typ_pos = -1;
2006  wvhdl_pos = -1;
2007
2008  if (r->typ==NULL)
2009    return FALSE;
2010
2011
2012  int j = p; // Which IS record to use...
2013  for( int pos = 0; pos < r->OrdSize; pos++ )
2014    if( r->typ[pos].ord_typ == ro_am)
2015      if( j-- == 0 )
2016      {
2017        typ_pos = pos;
2018
2019        if( bSearchWvhdl )
2020        {
2021          const int nblocks = rBlocks(r) - 1;
2022          const int* w = r->typ[pos].data.am.weights; // ?
2023
2024          for( pos = 0; pos <= nblocks; pos ++ )
2025            if (r->order[pos] == ringorder_am)
2026              if( r->wvhdl[pos] == w )
2027              {
2028                wvhdl_pos = pos;
2029                break;
2030              }
2031          if (wvhdl_pos < 0)
2032            return FALSE;
2033
2034          assume(wvhdl_pos >= 0);
2035        }
2036        assume(typ_pos >= 0);
2037        return TRUE;
2038      }
2039
2040  return FALSE;
2041}
2042
2043// // ?
2044// static BOOLEAN GetAMData(leftv res, leftv h)
2045// {
2046//   NoReturn(res);
2047//
2048//   const ring r = currRing;
2049//
2050//   int p = 0; // which IS-block? p^th!
2051//
2052//   if ((h!=NULL) && (h->Typ()==INT_CMD))
2053//     p = (int)((long)(h->Data())); h=h->next;
2054//
2055//   assume(p >= 0);
2056//
2057//   int d, w;
2058//   
2059//   if( !rGetAMPos(r, p, d, w, TRUE) )
2060//   {
2061//     Werror("`GetAMData([int])`: no %d^th _am block-ordering!", p);
2062//     return TRUE;
2063//   }
2064//
2065//   assume( r->typ[d].ord_typ == ro_am );
2066//   assume( r->order[w] == ringorder_am );
2067//
2068//
2069//   const short start = r->typ[d].data.am.start;  // bounds of ordering (in E)
2070//   const short end = r->typ[d].data.am.end;
2071//   const short len_gen = r->typ[d].data.am.len_gen; // i>len_gen: weight(gen(i)):=0
2072//   const int *weights = r->typ[d].data.am.weights; // pointers into wvhdl field of length (end-start+1) + len_gen
2073//   // contents w_1,... w_n, len, mod_w_1, .. mod_w_len, 0
2074//
2075//   assume( weights == r->wvhdl[w] );
2076//
2077//   
2078//   lists l=(lists)omAllocBin(slists_bin);
2079//   l->Init(2);
2080//
2081//   const short V = end-start+1;
2082//   intvec* ww_vars = new intvec(V);
2083//   intvec* ww_gens = new intvec(len_gen);
2084//
2085//   for (int i = 0; i < V; i++ )
2086//     (*ww_vars)[i] = weights[i];
2087//
2088//   assume( weights[V] == len_gen );
2089//
2090//   for (int i = 0; i < len_gen; i++ )
2091//     (*ww_gens)[i] = weights[i - V - 1];
2092//   
2093//
2094//   l->m[0].rtyp = INTVEC_CMD;
2095//   l->m[0].data = reinterpret_cast<void *>(ww_vars);
2096//
2097//   l->m[1].rtyp = INTVEC_CMD;
2098//   l->m[1].data = reinterpret_cast<void *>(ww_gens);
2099//
2100//
2101//   return FALSE;
2102//
2103// }
2104*/
2105
2106/// Returns old SyzCompLimit, can set new limit
2107static BOOLEAN SetInducedReferrence(leftv res, leftv h)
2108{
2109  NoReturn(res);
2110
2111  const ring r = currRing;
2112
2113  if( !( (h!=NULL) && ( (h->Typ()==IDEAL_CMD) || (h->Typ()==MODUL_CMD))) )
2114  {
2115    WerrorS("`SetInducedReferrence(<ideal/module>, [int[, int]])` expected");
2116    return TRUE;
2117  }
2118
2119  const ideal F = (ideal)h->Data(); ; // No copy!
2120  h=h->next;
2121
2122  int rank = 0;
2123
2124  if ((h!=NULL) && (h->Typ()==INT_CMD))
2125  {
2126    rank = (int)((long)(h->Data())); h=h->next;
2127    assume(rank >= 0);
2128  } else
2129    rank = id_RankFreeModule(F, r); // Starting syz-comp (1st: i+1)
2130
2131  int p = 0; // which IS-block? p^th!
2132
2133  if ((h!=NULL) && (h->Typ()==INT_CMD))
2134  {
2135    p = (int)((long)(h->Data())); h=h->next;
2136    assume(p >= 0);
2137  }
2138
2139  const int posIS = rGetISPos(p, r);
2140
2141  if(  /*(*/ -1 == posIS /*)*/  )
2142  {
2143    WerrorS("`SetInducedReferrence(<ideal/module>, [int[, int]])` called on incompatible ring (not created by 'MakeInducedSchreyerOrdering'!)");
2144    return TRUE;
2145  }
2146
2147
2148
2149  // F & componentWeights belong to that ordering block of currRing now:
2150  rSetISReference(r, F, rank, p); // F will be copied!
2151  return FALSE;
2152}
2153
2154
2155//    F = ISUpdateComponents( F, V, MIN );
2156//    // replace gen(i) -> gen(MIN + V[i-MIN]) for all i > MIN in all terms from F!
2157static BOOLEAN ISUpdateComponents(leftv res, leftv h)
2158{
2159  NoReturn(res);
2160
2161  PrintS("ISUpdateComponents:.... \n");
2162
2163  if ((h!=NULL) && (h->Typ()==MODUL_CMD))
2164  {
2165    ideal F = (ideal)h->Data(); ; // No copy!
2166    h=h->next;
2167
2168    if ((h!=NULL) && (h->Typ()==INTVEC_CMD))
2169    {
2170      const intvec* const V = (const intvec* const) h->Data();
2171      h=h->next;
2172
2173      if ((h!=NULL) && (h->Typ()==INT_CMD))
2174      {
2175        const int MIN = (int)((long)(h->Data()));
2176
2177        pISUpdateComponents(F, V, MIN, currRing);
2178        return FALSE;
2179      }
2180    }
2181  }
2182
2183  WerrorS("`ISUpdateComponents(<module>, intvec, int)` expected");
2184  return TRUE;
2185}
2186
2187
2188/// NF using length
2189static BOOLEAN reduce_syz(leftv res, leftv h)
2190{
2191  // const ring r = currRing;
2192
2193  if ( !( (h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD) ) )
2194  {
2195    WerrorS("`reduce_syz(<poly/vector>!, <ideal/module>, <int>, [int])` expected");
2196    return TRUE;
2197  }
2198
2199  res->rtyp = h->Typ();
2200  const poly v = reinterpret_cast<poly>(h->Data());
2201  h=h->next;
2202
2203  if ( !( (h!=NULL) && (h->Typ()==MODUL_CMD || h->Typ()==IDEAL_CMD ) ) )
2204  {
2205    WerrorS("`reduce_syz(<poly/vector>, <ideal/module>!, <int>, [int])` expected");
2206    return TRUE;
2207  }
2208
2209  assumeStdFlag(h);
2210  const ideal M = reinterpret_cast<ideal>(h->Data()); h=h->next;
2211
2212
2213  if ( !( (h!=NULL) && (h->Typ()== INT_CMD)  ) )
2214  {
2215    WerrorS("`reduce_syz(<poly/vector>, <ideal/module>, <int>!, [int])` expected");
2216    return TRUE;
2217  }
2218
2219  const int iSyzComp = (int)((long)(h->Data())); h=h->next;
2220
2221  int iLazyReduce = 0;
2222
2223  if ( ( (h!=NULL) && (h->Typ()== INT_CMD)  ) )
2224    iLazyReduce = (int)((long)(h->Data())); 
2225
2226  res->data = (void *)kNFLength(M, currQuotient, v, iSyzComp, iLazyReduce); // NOTE: currRing :(
2227  return FALSE;
2228}
2229
2230
2231/// Get raw syzygies (idPrepare)
2232static BOOLEAN idPrepare(leftv res, leftv h)
2233{
2234  //        extern int rGetISPos(const int p, const ring r);
2235
2236  const ring r = currRing;
2237
2238  const bool isSyz = rIsSyzIndexRing(r);
2239  const int posIS = rGetISPos(0, r);
2240
2241
2242  if ( !( (h!=NULL) && (h->Typ()==MODUL_CMD) && (h->Data() != NULL) ) )
2243  {
2244    WerrorS("`idPrepare(<module>)` expected");
2245    return TRUE;
2246  }
2247
2248  const ideal I = reinterpret_cast<ideal>(h->Data());
2249
2250  assume( I != NULL );
2251  idTest(I);
2252
2253  int iComp = -1;
2254
2255  h=h->next;
2256  if ( (h!=NULL) && (h->Typ()==INT_CMD) )
2257  {
2258    iComp = (int)((long)(h->Data()));
2259  } else
2260  {
2261      if( (!isSyz) && (-1 == posIS) )
2262      {
2263        WerrorS("`idPrepare(<...>)` called on incompatible ring (not created by 'MakeSyzCompOrdering' or 'MakeInducedSchreyerOrdering'!)");
2264        return TRUE;
2265      }
2266
2267    if( isSyz )
2268      iComp = rGetCurrSyzLimit(r);
2269    else
2270      iComp = id_RankFreeModule(r->typ[posIS].data.is.F, r); // ;
2271  }
2272 
2273  assume(iComp >= 0);
2274
2275
2276  intvec* w = reinterpret_cast<intvec *>(atGet(h, "isHomog", INTVEC_CMD));
2277  tHomog hom = testHomog;
2278
2279  //           int add_row_shift = 0;
2280  //
2281  if (w!=NULL)
2282  {
2283    w = ivCopy(w);
2284  //             add_row_shift = ww->min_in();
2285  //
2286  //             (*ww) -= add_row_shift;
2287  //             
2288  //             if (idTestHomModule(I, currQuotient, ww))
2289  //             {
2290    hom = isHomog;
2291  //               w = ww;
2292  //             }
2293  //             else
2294  //             {
2295  //               //WarnS("wrong weights");
2296  //               delete ww;
2297  //               w = NULL;
2298  //               hom=testHomog;
2299  //             }
2300  }
2301
2302
2303  // computes syzygies of h1,
2304  // works always in a ring with ringorder_s
2305  // NOTE: rSetSyzComp(syzcomp) should better be called beforehand
2306  //        ideal idPrepare (ideal  h1, tHomog hom, int syzcomp, intvec **w);
2307
2308  ideal J = // idPrepare( I, hom, iComp, &w);
2309           kStd(I, currQuotient, hom, &w, NULL, iComp);
2310
2311  idTest(J);
2312
2313  if (w!=NULL)
2314    atSet(res, omStrDup("isHomog"), w, INTVEC_CMD);
2315  //             if (w!=NULL) delete w;
2316
2317  res->rtyp = MODUL_CMD;
2318  res->data = reinterpret_cast<void *>(J);
2319  return FALSE;
2320}
2321
2322/// Get raw syzygies (idPrepare)
2323static BOOLEAN _p_Content(leftv res, leftv h)
2324{
2325  if ( !( (h!=NULL) && (h->Typ()==POLY_CMD) && (h->Data() != NULL) ) )
2326  {
2327    WerrorS("`p_Content(<poly-var>)` expected");
2328    return TRUE;
2329  }
2330
2331
2332  const poly p = reinterpret_cast<poly>(h->Data());
2333
2334 
2335  pTest(p);  pWrite(p); PrintLn();
2336
2337 
2338  p_Content( p, currRing);     
2339
2340  pTest(p);
2341  pWrite(p); PrintLn();
2342 
2343  NoReturn(res);
2344  return FALSE;
2345}
2346
2347static BOOLEAN _m2_end(leftv res, leftv h)
2348{
2349  int ret = 0;
2350 
2351  if ( (h!=NULL) && (h->Typ()!=INT_CMD) )
2352  {
2353    WerrorS("`m2_end([<int>])` expected");
2354    return TRUE;
2355  }
2356  ret = (int)(long)(h->Data());
2357
2358  m2_end( ret );
2359
2360  NoReturn(res);
2361  return FALSE;
2362}
2363
2364   
2365
2366END_NAMESPACE
2367
2368
2369int SI_MOD_INIT(syzextra)(SModulFunctions* psModulFunctions) 
2370{
2371#define ADD0(A,B,C,D,E) A(B, (char*)C, D, E)
2372// #define ADD(A,B,C,D,E) ADD0(iiAddCproc, "", C, D, E)
2373  #define ADD(A,B,C,D,E) ADD0(A->iiAddCproc, B, C, D, E)
2374  ADD(psModulFunctions, currPack->libname, "ClearContent", FALSE, _ClearContent);
2375  ADD(psModulFunctions, currPack->libname, "ClearDenominators", FALSE, _ClearDenominators);
2376
2377  ADD(psModulFunctions, currPack->libname, "m2_end", FALSE, _m2_end);
2378
2379  ADD(psModulFunctions, currPack->libname, "DetailedPrint", FALSE, DetailedPrint);
2380  ADD(psModulFunctions, currPack->libname, "leadmonomial", FALSE, leadmonom);
2381  ADD(psModulFunctions, currPack->libname, "leadcomp", FALSE, leadcomp);
2382  ADD(psModulFunctions, currPack->libname, "leadrawexp", FALSE, leadrawexp);
2383
2384  ADD(psModulFunctions, currPack->libname, "ISUpdateComponents", FALSE, ISUpdateComponents);
2385  ADD(psModulFunctions, currPack->libname, "SetInducedReferrence", FALSE, SetInducedReferrence);
2386  ADD(psModulFunctions, currPack->libname, "GetInducedData", FALSE, GetInducedData);
2387  ADD(psModulFunctions, currPack->libname, "SetSyzComp", FALSE, SetSyzComp);
2388  ADD(psModulFunctions, currPack->libname, "MakeInducedSchreyerOrdering", FALSE, MakeInducedSchreyerOrdering);
2389  ADD(psModulFunctions, currPack->libname, "MakeSyzCompOrdering", FALSE, MakeSyzCompOrdering);
2390
2391  ADD(psModulFunctions, currPack->libname, "ProfilerStart", FALSE, _ProfilerStart); ADD(psModulFunctions, currPack->libname, "ProfilerStop",  FALSE, _ProfilerStop );
2392 
2393  ADD(psModulFunctions, currPack->libname, "noop", FALSE, noop);
2394  ADD(psModulFunctions, currPack->libname, "idPrepare", FALSE, idPrepare);
2395  ADD(psModulFunctions, currPack->libname, "reduce_syz", FALSE, reduce_syz);
2396
2397  ADD(psModulFunctions, currPack->libname, "p_Content", FALSE, _p_Content);
2398
2399  ADD(psModulFunctions, currPack->libname, "Tail", FALSE, Tail);
2400 
2401  ADD(psModulFunctions, currPack->libname, "ComputeLeadingSyzygyTerms", FALSE, ComputeLeadingSyzygyTerms);
2402  ADD(psModulFunctions, currPack->libname, "Compute2LeadingSyzygyTerms", FALSE, Compute2LeadingSyzygyTerms);
2403 
2404  ADD(psModulFunctions, currPack->libname, "Sort_c_ds", FALSE, Sort_c_ds);
2405  ADD(psModulFunctions, currPack->libname, "FindReducer", FALSE, FindReducer);
2406
2407
2408  ADD(psModulFunctions, currPack->libname, "ReduceTerm", FALSE, ReduceTerm);
2409  ADD(psModulFunctions, currPack->libname, "TraverseTail", FALSE, TraverseTail);
2410
2411   
2412  ADD(psModulFunctions, currPack->libname, "SchreyerSyzygyNF", FALSE, SchreyerSyzygyNF);
2413 
2414  //  ADD(psModulFunctions, currPack->libname, "GetAMData", FALSE, GetAMData);
2415
2416  //  ADD(psModulFunctions, currPack->libname, "", FALSE, );
2417
2418#undef ADD 
2419  return 0;
2420}
2421
2422#ifndef EMBED_PYTHON
2423extern "C" { 
2424int mod_init(SModulFunctions* psModulFunctions)
2425{ 
2426  return SI_MOD_INIT(syzextra)(psModulFunctions); 
2427}
2428}
2429#endif
Note: See TracBrowser for help on using the repository browser.