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

spielwiese
Last change on this file since f6c459 was f6c459, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
Cleanup and minor changes
  • Property mode set to 100644
File size: 44.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
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/// TODO: save shortcut (syz: |-.->) LM(LM(m) * "t") -> syz?
1100/// proc SSFindReducer(def product, def syzterm, def L, def T, list #)
1101static BOOLEAN FindReducer(leftv res, leftv h)
1102{
1103  const char* usage = "`FindReducer(<poly/vector>, <vector/0>, <ideal/module>[,<module>])` expected";
1104  const ring r = currRing;
1105
1106  NoReturn(res);
1107
1108
1109#ifndef NDEBUG
1110  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
1111#else
1112  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
1113#endif
1114
1115  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
1116
1117  if ((h==NULL) || (h->Typ()!=VECTOR_CMD && h->Typ() !=POLY_CMD) || (h->Data() == NULL))
1118  {
1119    WerrorS(usage);
1120    return TRUE;   
1121  }
1122   
1123  const poly product = (poly) h->Data(); assume (product != NULL);
1124
1125
1126  h = h->Next();
1127  if ((h==NULL) || !((h->Typ()==VECTOR_CMD) || (h->Data() == NULL)) )
1128  {
1129    WerrorS(usage);
1130    return TRUE;   
1131  }
1132
1133  const poly syzterm = (poly) h->Data(); h = h->Next();
1134
1135  int c = 0;
1136 
1137  if (syzterm != NULL)
1138    c = p_GetComp(syzterm, r) - 1;
1139 
1140 
1141  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1142  {
1143    WerrorS(usage);
1144    return TRUE;   
1145  }
1146 
1147  const ideal L = (ideal) h->Data(); h = h->Next();
1148
1149  assume( IDELEMS(L) > 0 );
1150  assume( c >= 0 && c < IDELEMS(L) );
1151
1152 
1153/*
1154  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
1155  {
1156    WerrorS(usage);
1157    return TRUE;   
1158  }
1159
1160  const ideal T = (ideal) h->Data(); h = h->Next();
1161*/
1162
1163  ideal LS = NULL;
1164
1165  if ((h != NULL) && (h->Typ() ==MODUL_CMD) && (h->Data() != NULL))
1166  {
1167    LS = (ideal)h->Data();
1168    h = h->Next();
1169  }
1170
1171  if( __TAILREDSYZ__ )
1172    assume (LS != NULL);
1173
1174  assume( h == NULL );
1175
1176  if( __DEBUG__ )
1177  {
1178    PrintS("FindReducer(product, syzterm, L, T, #)::Input: \n");
1179
1180    PrintS("product: "); dPrint(product, r, r, 2);
1181    PrintS("syzterm: "); dPrint(syzterm, r, r, 2);
1182    PrintS("L: "); dPrint(L, r, r, 2);
1183//    PrintS("T: "); dPrint(T, r, r, 4);
1184
1185    if( LS == NULL )
1186      PrintS("LS: NULL\n");
1187    else
1188    {
1189      PrintS("LS: "); dPrint(LS, r, r, 2);
1190    }
1191  }
1192
1193
1194  if (__SYZCHECK__ && syzterm != NULL)
1195  {
1196    const poly m = L->m[c];
1197
1198    assume( m != NULL ); assume( pNext(m) == NULL );
1199
1200    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
1201    assume( p_EqualPolys(lm, product, r) );
1202
1203    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
1204    //  def @@product = leadmonomial(syzterm) * L[@@r];
1205
1206    p_Delete(&lm, r);   
1207  }
1208
1209 
1210  res->rtyp = VECTOR_CMD;
1211
1212  // looking for an appropriate diviser q = L[k]...
1213  for( int k = IDELEMS(L)-1; k>= 0; k-- )
1214  {
1215    const poly p = L->m[k];   
1216
1217    // ... which divides the product, looking for the _1st_ appropriate one!
1218    if( !p_LmDivisibleBy(p, product, r) )
1219      continue;
1220
1221
1222    const poly q = p_New(r);
1223    pNext(q) = NULL;
1224    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
1225    p_SetComp(q, k + 1, r);
1226    p_Setm(q, r);
1227
1228    // cannot allow something like: a*gen(i) - a*gen(i)
1229    if (syzterm != NULL && (k == c))
1230      if (p_ExpVectorEqual(syzterm, q, r))
1231      {
1232        if( __DEBUG__ )
1233        {
1234          Print("FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1235          dPrint(syzterm, r, r, 1);
1236        }   
1237       
1238        p_LmFree(q, r);
1239        continue;
1240      }
1241
1242    // while the complement (the fraction) is not reducible by leading syzygies
1243    if( LS != NULL )
1244    {
1245      BOOLEAN ok = TRUE;
1246     
1247      for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )
1248      {
1249        const poly pp = LS->m[kk];
1250
1251        if( p_LmDivisibleBy(pp, q, r) )
1252        {
1253          if( __DEBUG__ )
1254          {
1255            Print("FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);
1256            dPrint(pp, r, r, 1);
1257          }   
1258
1259          ok = FALSE; // q in <LS> :((
1260          break;                 
1261        }
1262      }
1263     
1264      if(!ok)
1265      {
1266        p_LmFree(q, r);
1267        continue;
1268      }
1269    }
1270
1271    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
1272
1273    if( __DEBUG__ )
1274    {
1275      PrintS("FindReducer::Output: \n");
1276      dPrint(q, r, r, 1);
1277    }   
1278
1279    res->data = q;
1280    return FALSE;
1281  }
1282
1283  res->data = NULL;
1284  return FALSE;   
1285 
1286}
1287
1288
1289
1290/// Get leading term without a module component
1291static BOOLEAN leadmonom(leftv res, leftv h)
1292{
1293  NoReturn(res);
1294
1295  if ((h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD) && (h->Data() != NULL))
1296  {
1297    const ring r = currRing;
1298    const poly p = (poly)(h->Data());
1299
1300    const poly m = leadmonom(p, r);
1301
1302    res->rtyp = POLY_CMD;
1303    res->data = reinterpret_cast<void *>(m);
1304
1305    return FALSE;
1306  }
1307
1308  WerrorS("`leadmonom(<poly/vector>)` expected");
1309  return TRUE;
1310}
1311
1312/// Get leading component
1313static BOOLEAN leadcomp(leftv res, leftv h)
1314{
1315  NoReturn(res);
1316
1317  if ((h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD))
1318  {
1319    const ring r = currRing;
1320
1321    const poly p = (poly)(h->Data());
1322
1323    if (p != NULL )
1324    {
1325      assume( p != NULL );
1326      assume( p_LmTest(p, r) );
1327
1328      const unsigned long iComp = p_GetComp(p, r);
1329
1330  //    assume( iComp > 0 ); // p is a vector
1331
1332      res->data = reinterpret_cast<void *>(jjLONG2N(iComp));
1333    } else
1334      res->data = reinterpret_cast<void *>(jjLONG2N(0));
1335     
1336
1337    res->rtyp = BIGINT_CMD;
1338    return FALSE;
1339  }
1340
1341  WerrorS("`leadcomp(<poly/vector>)` expected");
1342  return TRUE;
1343}
1344
1345
1346
1347
1348/// Get raw leading exponent vector
1349static BOOLEAN leadrawexp(leftv res, leftv h)
1350{
1351  NoReturn(res);
1352
1353  if ((h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD) && (h->Data() != NULL))
1354  {
1355    const ring r = currRing;
1356    const poly p = (poly)(h->Data());
1357
1358    assume( p != NULL );
1359    assume( p_LmTest(p, r) );
1360
1361    const int iExpSize = r->ExpL_Size;
1362
1363//    intvec *iv = new intvec(iExpSize);
1364
1365    lists l=(lists)omAllocBin(slists_bin);
1366    l->Init(iExpSize);
1367
1368    for(int i = iExpSize-1; i >= 0; i--)
1369    {
1370      l->m[i].rtyp = BIGINT_CMD;
1371      l->m[i].data = reinterpret_cast<void *>(jjLONG2N(p->exp[i])); // longs...
1372    }
1373
1374    res->rtyp = LIST_CMD; // list of bigints
1375    res->data = reinterpret_cast<void *>(l);
1376    return FALSE;
1377  }
1378
1379  WerrorS("`leadrawexp(<poly/vector>)` expected");
1380  return TRUE;
1381}
1382
1383
1384/// Endowe the current ring with additional (leading) Syz-component ordering
1385static BOOLEAN MakeSyzCompOrdering(leftv res, leftv /*h*/)
1386{
1387
1388  NoReturn(res);
1389
1390  //    res->data = rCurrRingAssure_SyzComp(); // changes current ring! :(
1391  res->data = reinterpret_cast<void *>(rAssure_SyzComp(currRing, TRUE));
1392  res->rtyp = RING_CMD; // return new ring!
1393  // QRING_CMD?
1394
1395  return FALSE;
1396}
1397
1398
1399/// Same for Induced Schreyer ordering (ordering on components is defined by sign!)
1400static BOOLEAN MakeInducedSchreyerOrdering(leftv res, leftv h)
1401{
1402
1403  NoReturn(res);
1404
1405  int sign = 1;
1406  if ((h!=NULL) && (h->Typ()==INT_CMD))
1407  {
1408    const int s = (int)((long)(h->Data()));
1409
1410    if( s != -1 && s != 1 )
1411    {
1412      WerrorS("`MakeInducedSchreyerOrdering(<int>)` called with wrong integer argument (must be +-1)!");
1413      return TRUE;
1414    }
1415
1416    sign = s;           
1417  }
1418
1419  assume( sign == 1 || sign == -1 );
1420  res->data = reinterpret_cast<void *>(rAssure_InducedSchreyerOrdering(currRing, TRUE, sign));
1421  res->rtyp = RING_CMD; // return new ring!
1422  // QRING_CMD?
1423  return FALSE;
1424}
1425
1426
1427/// Returns old SyzCompLimit, can set new limit
1428static BOOLEAN SetSyzComp(leftv res, leftv h)
1429{
1430  NoReturn(res);
1431
1432  const ring r = currRing;
1433
1434  if( !rIsSyzIndexRing(r) )
1435  {
1436    WerrorS("`SetSyzComp(<int>)` called on incompatible ring (not created by 'MakeSyzCompOrdering'!)");
1437    return TRUE;
1438  }
1439
1440  res->rtyp = INT_CMD;
1441  res->data = reinterpret_cast<void *>(rGetCurrSyzLimit(r)); // return old syz limit
1442
1443  if ((h!=NULL) && (h->Typ()==INT_CMD))
1444  {
1445    const int iSyzComp = (int)reinterpret_cast<long>(h->Data());
1446    assume( iSyzComp > 0 );
1447    rSetSyzComp(iSyzComp, currRing);
1448  }
1449
1450  return FALSE;
1451}
1452
1453/// ?
1454static BOOLEAN GetInducedData(leftv res, leftv h)
1455{
1456  NoReturn(res);
1457
1458  const ring r = currRing;
1459
1460  int p = 0; // which IS-block? p^th!
1461
1462  if ((h!=NULL) && (h->Typ()==INT_CMD))
1463  {
1464    p = (int)((long)(h->Data())); h=h->next;
1465    assume(p >= 0);
1466  }
1467
1468  const int pos = rGetISPos(p, r);
1469
1470  if(  /*(*/ -1 == pos /*)*/  )
1471  {
1472    WerrorS("`GetInducedData([int])` called on incompatible ring (not created by 'MakeInducedSchreyerOrdering'!)");
1473    return TRUE;
1474  }
1475
1476
1477  const int iLimit = r->typ[pos].data.is.limit;
1478  const ideal F = r->typ[pos].data.is.F;
1479  ideal FF = id_Copy(F, r);
1480
1481  lists l=(lists)omAllocBin(slists_bin);
1482  l->Init(2);
1483
1484  l->m[0].rtyp = INT_CMD;
1485  l->m[0].data = reinterpret_cast<void *>(iLimit);
1486
1487
1488  //        l->m[1].rtyp = MODUL_CMD;
1489
1490  if( idIsModule(FF, r) )
1491  {
1492    l->m[1].rtyp = MODUL_CMD;
1493
1494    //          Print("before: %d\n", FF->nrows);
1495    //          FF->nrows = id_RankFreeModule(FF, r); // ???
1496    //          Print("after: %d\n", FF->nrows);
1497  }
1498  else
1499    l->m[1].rtyp = IDEAL_CMD;
1500
1501  l->m[1].data = reinterpret_cast<void *>(FF);
1502
1503  res->rtyp = LIST_CMD; // list of int/module
1504  res->data = reinterpret_cast<void *>(l);
1505
1506  return FALSE;
1507
1508}
1509
1510
1511/* // the following turned out to be unnecessary...   
1512/// Finds p^th AM ordering, and returns its position in r->typ[] AND
1513/// corresponding &r->wvhdl[]
1514/// returns FALSE if something went wrong!
1515/// p - starts with 0!
1516BOOLEAN rGetAMPos(const ring r, const int p, int &typ_pos, int &wvhdl_pos, const BOOLEAN bSearchWvhdl = FALSE)
1517{
1518#if MYTEST
1519  Print("rGetAMPos(p: %d...)\nF:", p);
1520  PrintLn();
1521#endif
1522  typ_pos = -1;
1523  wvhdl_pos = -1;
1524
1525  if (r->typ==NULL)
1526    return FALSE;
1527
1528
1529  int j = p; // Which IS record to use...
1530  for( int pos = 0; pos < r->OrdSize; pos++ )
1531    if( r->typ[pos].ord_typ == ro_am)
1532      if( j-- == 0 )
1533      {
1534        typ_pos = pos;
1535
1536        if( bSearchWvhdl )
1537        {
1538          const int nblocks = rBlocks(r) - 1;
1539          const int* w = r->typ[pos].data.am.weights; // ?
1540
1541          for( pos = 0; pos <= nblocks; pos ++ )
1542            if (r->order[pos] == ringorder_am)
1543              if( r->wvhdl[pos] == w )
1544              {
1545                wvhdl_pos = pos;
1546                break;
1547              }
1548          if (wvhdl_pos < 0)
1549            return FALSE;
1550
1551          assume(wvhdl_pos >= 0);
1552        }
1553        assume(typ_pos >= 0);
1554        return TRUE;
1555      }
1556
1557  return FALSE;
1558}
1559
1560// // ?
1561// static BOOLEAN GetAMData(leftv res, leftv h)
1562// {
1563//   NoReturn(res);
1564//
1565//   const ring r = currRing;
1566//
1567//   int p = 0; // which IS-block? p^th!
1568//
1569//   if ((h!=NULL) && (h->Typ()==INT_CMD))
1570//     p = (int)((long)(h->Data())); h=h->next;
1571//
1572//   assume(p >= 0);
1573//
1574//   int d, w;
1575//   
1576//   if( !rGetAMPos(r, p, d, w, TRUE) )
1577//   {
1578//     Werror("`GetAMData([int])`: no %d^th _am block-ordering!", p);
1579//     return TRUE;
1580//   }
1581//
1582//   assume( r->typ[d].ord_typ == ro_am );
1583//   assume( r->order[w] == ringorder_am );
1584//
1585//
1586//   const short start = r->typ[d].data.am.start;  // bounds of ordering (in E)
1587//   const short end = r->typ[d].data.am.end;
1588//   const short len_gen = r->typ[d].data.am.len_gen; // i>len_gen: weight(gen(i)):=0
1589//   const int *weights = r->typ[d].data.am.weights; // pointers into wvhdl field of length (end-start+1) + len_gen
1590//   // contents w_1,... w_n, len, mod_w_1, .. mod_w_len, 0
1591//
1592//   assume( weights == r->wvhdl[w] );
1593//
1594//   
1595//   lists l=(lists)omAllocBin(slists_bin);
1596//   l->Init(2);
1597//
1598//   const short V = end-start+1;
1599//   intvec* ww_vars = new intvec(V);
1600//   intvec* ww_gens = new intvec(len_gen);
1601//
1602//   for (int i = 0; i < V; i++ )
1603//     (*ww_vars)[i] = weights[i];
1604//
1605//   assume( weights[V] == len_gen );
1606//
1607//   for (int i = 0; i < len_gen; i++ )
1608//     (*ww_gens)[i] = weights[i - V - 1];
1609//   
1610//
1611//   l->m[0].rtyp = INTVEC_CMD;
1612//   l->m[0].data = reinterpret_cast<void *>(ww_vars);
1613//
1614//   l->m[1].rtyp = INTVEC_CMD;
1615//   l->m[1].data = reinterpret_cast<void *>(ww_gens);
1616//
1617//
1618//   return FALSE;
1619//
1620// }
1621*/
1622
1623/// Returns old SyzCompLimit, can set new limit
1624static BOOLEAN SetInducedReferrence(leftv res, leftv h)
1625{
1626  NoReturn(res);
1627
1628  const ring r = currRing;
1629
1630  if( !( (h!=NULL) && ( (h->Typ()==IDEAL_CMD) || (h->Typ()==MODUL_CMD))) )
1631  {
1632    WerrorS("`SetInducedReferrence(<ideal/module>, [int[, int]])` expected");
1633    return TRUE;
1634  }
1635
1636  const ideal F = (ideal)h->Data(); ; // No copy!
1637  h=h->next;
1638
1639  int rank = 0;
1640
1641  if ((h!=NULL) && (h->Typ()==INT_CMD))
1642  {
1643    rank = (int)((long)(h->Data())); h=h->next;
1644    assume(rank >= 0);
1645  } else
1646    rank = id_RankFreeModule(F, r); // Starting syz-comp (1st: i+1)
1647
1648  int p = 0; // which IS-block? p^th!
1649
1650  if ((h!=NULL) && (h->Typ()==INT_CMD))
1651  {
1652    p = (int)((long)(h->Data())); h=h->next;
1653    assume(p >= 0);
1654  }
1655
1656  const int posIS = rGetISPos(p, r);
1657
1658  if(  /*(*/ -1 == posIS /*)*/  )
1659  {
1660    WerrorS("`SetInducedReferrence(<ideal/module>, [int[, int]])` called on incompatible ring (not created by 'MakeInducedSchreyerOrdering'!)");
1661    return TRUE;
1662  }
1663
1664
1665
1666  // F & componentWeights belong to that ordering block of currRing now:
1667  rSetISReference(r, F, rank, p); // F will be copied!
1668  return FALSE;
1669}
1670
1671
1672//    F = ISUpdateComponents( F, V, MIN );
1673//    // replace gen(i) -> gen(MIN + V[i-MIN]) for all i > MIN in all terms from F!
1674static BOOLEAN ISUpdateComponents(leftv res, leftv h)
1675{
1676  NoReturn(res);
1677
1678  PrintS("ISUpdateComponents:.... \n");
1679
1680  if ((h!=NULL) && (h->Typ()==MODUL_CMD))
1681  {
1682    ideal F = (ideal)h->Data(); ; // No copy!
1683    h=h->next;
1684
1685    if ((h!=NULL) && (h->Typ()==INTVEC_CMD))
1686    {
1687      const intvec* const V = (const intvec* const) h->Data();
1688      h=h->next;
1689
1690      if ((h!=NULL) && (h->Typ()==INT_CMD))
1691      {
1692        const int MIN = (int)((long)(h->Data()));
1693
1694        pISUpdateComponents(F, V, MIN, currRing);
1695        return FALSE;
1696      }
1697    }
1698  }
1699
1700  WerrorS("`ISUpdateComponents(<module>, intvec, int)` expected");
1701  return TRUE;
1702}
1703
1704
1705/// NF using length
1706static BOOLEAN reduce_syz(leftv res, leftv h)
1707{
1708  // const ring r = currRing;
1709
1710  if ( !( (h!=NULL) && (h->Typ()==VECTOR_CMD || h->Typ()==POLY_CMD) ) )
1711  {
1712    WerrorS("`reduce_syz(<poly/vector>!, <ideal/module>, <int>, [int])` expected");
1713    return TRUE;
1714  }
1715
1716  res->rtyp = h->Typ();
1717  const poly v = reinterpret_cast<poly>(h->Data());
1718  h=h->next;
1719
1720  if ( !( (h!=NULL) && (h->Typ()==MODUL_CMD || h->Typ()==IDEAL_CMD ) ) )
1721  {
1722    WerrorS("`reduce_syz(<poly/vector>, <ideal/module>!, <int>, [int])` expected");
1723    return TRUE;
1724  }
1725
1726  assumeStdFlag(h);
1727  const ideal M = reinterpret_cast<ideal>(h->Data()); h=h->next;
1728
1729
1730  if ( !( (h!=NULL) && (h->Typ()== INT_CMD)  ) )
1731  {
1732    WerrorS("`reduce_syz(<poly/vector>, <ideal/module>, <int>!, [int])` expected");
1733    return TRUE;
1734  }
1735
1736  const int iSyzComp = (int)((long)(h->Data())); h=h->next;
1737
1738  int iLazyReduce = 0;
1739
1740  if ( ( (h!=NULL) && (h->Typ()== INT_CMD)  ) )
1741    iLazyReduce = (int)((long)(h->Data())); 
1742
1743  res->data = (void *)kNFLength(M, currQuotient, v, iSyzComp, iLazyReduce); // NOTE: currRing :(
1744  return FALSE;
1745}
1746
1747
1748/// Get raw syzygies (idPrepare)
1749static BOOLEAN idPrepare(leftv res, leftv h)
1750{
1751  //        extern int rGetISPos(const int p, const ring r);
1752
1753  const ring r = currRing;
1754
1755  const bool isSyz = rIsSyzIndexRing(r);
1756  const int posIS = rGetISPos(0, r);
1757
1758
1759  if ( !( (h!=NULL) && (h->Typ()==MODUL_CMD) && (h->Data() != NULL) ) )
1760  {
1761    WerrorS("`idPrepare(<module>)` expected");
1762    return TRUE;
1763  }
1764
1765  const ideal I = reinterpret_cast<ideal>(h->Data());
1766
1767  assume( I != NULL );
1768  idTest(I);
1769
1770  int iComp = -1;
1771
1772  h=h->next;
1773  if ( (h!=NULL) && (h->Typ()==INT_CMD) )
1774  {
1775    iComp = (int)((long)(h->Data()));
1776  } else
1777  {
1778      if( (!isSyz) && (-1 == posIS) )
1779      {
1780        WerrorS("`idPrepare(<...>)` called on incompatible ring (not created by 'MakeSyzCompOrdering' or 'MakeInducedSchreyerOrdering'!)");
1781        return TRUE;
1782      }
1783
1784    if( isSyz )
1785      iComp = rGetCurrSyzLimit(r);
1786    else
1787      iComp = id_RankFreeModule(r->typ[posIS].data.is.F, r); // ;
1788  }
1789 
1790  assume(iComp >= 0);
1791
1792
1793  intvec* w = reinterpret_cast<intvec *>(atGet(h, "isHomog", INTVEC_CMD));
1794  tHomog hom = testHomog;
1795
1796  //           int add_row_shift = 0;
1797  //
1798  if (w!=NULL)
1799  {
1800    w = ivCopy(w);
1801  //             add_row_shift = ww->min_in();
1802  //
1803  //             (*ww) -= add_row_shift;
1804  //             
1805  //             if (idTestHomModule(I, currQuotient, ww))
1806  //             {
1807    hom = isHomog;
1808  //               w = ww;
1809  //             }
1810  //             else
1811  //             {
1812  //               //WarnS("wrong weights");
1813  //               delete ww;
1814  //               w = NULL;
1815  //               hom=testHomog;
1816  //             }
1817  }
1818
1819
1820  // computes syzygies of h1,
1821  // works always in a ring with ringorder_s
1822  // NOTE: rSetSyzComp(syzcomp) should better be called beforehand
1823  //        ideal idPrepare (ideal  h1, tHomog hom, int syzcomp, intvec **w);
1824
1825  ideal J = // idPrepare( I, hom, iComp, &w);
1826           kStd(I, currQuotient, hom, &w, NULL, iComp);
1827
1828  idTest(J);
1829
1830  if (w!=NULL)
1831    atSet(res, omStrDup("isHomog"), w, INTVEC_CMD);
1832  //             if (w!=NULL) delete w;
1833
1834  res->rtyp = MODUL_CMD;
1835  res->data = reinterpret_cast<void *>(J);
1836  return FALSE;
1837}
1838
1839/// Get raw syzygies (idPrepare)
1840static BOOLEAN _p_Content(leftv res, leftv h)
1841{
1842  if ( !( (h!=NULL) && (h->Typ()==POLY_CMD) && (h->Data() != NULL) ) )
1843  {
1844    WerrorS("`p_Content(<poly-var>)` expected");
1845    return TRUE;
1846  }
1847
1848
1849  const poly p = reinterpret_cast<poly>(h->Data());
1850
1851 
1852  pTest(p);  pWrite(p); PrintLn();
1853
1854 
1855  p_Content( p, currRing);     
1856
1857  pTest(p);
1858  pWrite(p); PrintLn();
1859 
1860  NoReturn(res);
1861  return FALSE;
1862}
1863
1864static BOOLEAN _m2_end(leftv res, leftv h)
1865{
1866  int ret = 0;
1867 
1868  if ( (h!=NULL) && (h->Typ()!=INT_CMD) )
1869  {
1870    WerrorS("`m2_end([<int>])` expected");
1871    return TRUE;
1872  }
1873  ret = (int)(long)(h->Data());
1874
1875  m2_end( ret );
1876
1877  NoReturn(res);
1878  return FALSE;
1879}
1880
1881   
1882
1883END_NAMESPACE
1884
1885
1886int SI_MOD_INIT(syzextra)(SModulFunctions* psModulFunctions) 
1887{
1888#define ADD0(A,B,C,D,E) A(B, (char*)C, D, E)
1889// #define ADD(A,B,C,D,E) ADD0(iiAddCproc, "", C, D, E)
1890  #define ADD(A,B,C,D,E) ADD0(A->iiAddCproc, B, C, D, E)
1891  ADD(psModulFunctions, currPack->libname, "ClearContent", FALSE, _ClearContent);
1892  ADD(psModulFunctions, currPack->libname, "ClearDenominators", FALSE, _ClearDenominators);
1893
1894  ADD(psModulFunctions, currPack->libname, "m2_end", FALSE, _m2_end);
1895
1896  ADD(psModulFunctions, currPack->libname, "DetailedPrint", FALSE, DetailedPrint);
1897  ADD(psModulFunctions, currPack->libname, "leadmonomial", FALSE, leadmonom);
1898  ADD(psModulFunctions, currPack->libname, "leadcomp", FALSE, leadcomp);
1899  ADD(psModulFunctions, currPack->libname, "leadrawexp", FALSE, leadrawexp);
1900
1901  ADD(psModulFunctions, currPack->libname, "ISUpdateComponents", FALSE, ISUpdateComponents);
1902  ADD(psModulFunctions, currPack->libname, "SetInducedReferrence", FALSE, SetInducedReferrence);
1903  ADD(psModulFunctions, currPack->libname, "GetInducedData", FALSE, GetInducedData);
1904  ADD(psModulFunctions, currPack->libname, "SetSyzComp", FALSE, SetSyzComp);
1905  ADD(psModulFunctions, currPack->libname, "MakeInducedSchreyerOrdering", FALSE, MakeInducedSchreyerOrdering);
1906  ADD(psModulFunctions, currPack->libname, "MakeSyzCompOrdering", FALSE, MakeSyzCompOrdering);
1907
1908  ADD(psModulFunctions, currPack->libname, "ProfilerStart", FALSE, _ProfilerStart); ADD(psModulFunctions, currPack->libname, "ProfilerStop",  FALSE, _ProfilerStop );
1909 
1910  ADD(psModulFunctions, currPack->libname, "noop", FALSE, noop);
1911  ADD(psModulFunctions, currPack->libname, "idPrepare", FALSE, idPrepare);
1912  ADD(psModulFunctions, currPack->libname, "reduce_syz", FALSE, reduce_syz);
1913
1914  ADD(psModulFunctions, currPack->libname, "p_Content", FALSE, _p_Content);
1915
1916  ADD(psModulFunctions, currPack->libname, "Tail", FALSE, Tail);
1917 
1918  ADD(psModulFunctions, currPack->libname, "ComputeLeadingSyzygyTerms", FALSE, ComputeLeadingSyzygyTerms);
1919  ADD(psModulFunctions, currPack->libname, "Compute2LeadingSyzygyTerms", FALSE, Compute2LeadingSyzygyTerms);
1920 
1921  ADD(psModulFunctions, currPack->libname, "Sort_c_ds", FALSE, Sort_c_ds);
1922  ADD(psModulFunctions, currPack->libname, "FindReducer", FALSE, FindReducer);
1923 
1924 
1925  //  ADD(psModulFunctions, currPack->libname, "GetAMData", FALSE, GetAMData);
1926
1927  //  ADD(psModulFunctions, currPack->libname, "", FALSE, );
1928
1929#undef ADD 
1930  return 0;
1931}
1932
1933#ifndef EMBED_PYTHON
1934extern "C" { 
1935int mod_init(SModulFunctions* psModulFunctions)
1936{ 
1937  return SI_MOD_INIT(syzextra)(psModulFunctions); 
1938}
1939}
1940#endif
Note: See TracBrowser for help on using the repository browser.