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

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