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

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