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

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