source: git/kernel/nc.cc @ 5fe834

spielwiese
Last change on this file since 5fe834 was 5fe834, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: libpolys and kernel should use rDelete instead of rKill!
  • Property mode set to 100644
File size: 8.8 KB
Line 
1#include "mod2.h"
2#include <misc/auxiliary.h>
3
4#include <misc/options.h>
5
6#include <polys/simpleideals.h>
7#include <kernel/polys.h>
8
9#include <polys/prCopy.h>
10
11#include "ideals.h"
12#include "kstd1.h"
13
14#include "nc.h"
15
16ideal twostd(ideal I) // works in currRing only!
17{
18  ideal J = kStd(I, currQuotient, testHomog, NULL, NULL, 0, 0, NULL); // in currRing!!!
19  idSkipZeroes(J); // ring independent!
20
21  const int rN = currRing->N;
22
23  loop
24  {
25    ideal     K    = NULL;
26    const int s    = idElem(J); // ring independent
27
28    for(int i = 0; i < s; i++)
29    {
30      const poly p = J->m[i];
31
32#ifdef PDEBUG
33      p_Test(p, currRing);
34#if 0
35      Print("p: "); // !
36      p_Write(p, currRing);
37#endif
38#endif
39
40      for (int j = 1; j <= rN; j++) // for all j = 1..N
41      {
42        poly varj = p_One( currRing);
43        p_SetExp(varj, j, 1, currRing);
44        p_Setm(varj, currRing);
45
46        poly q = pp_Mult_mm(p, varj, currRing); // q = J[i] * var(j),
47
48#ifdef PDEBUG
49        p_Test(varj, currRing);
50        p_Test(p, currRing);
51        p_Test(q, currRing);
52#if 0
53        Print("Reducing p: "); // !
54        p_Write(p, currRing);
55        Print("With q: "); // !
56        p_Write(q, currRing);
57#endif
58#endif
59
60        p_Delete(&varj, currRing);
61
62        if (q != NULL)
63        {
64#ifdef PDEBUG
65#if 0
66          Print("Reducing q[j = %d]: ", j); // !
67          p_Write(q, currRing);
68
69          Print("With p:");
70          p_Write(p, currRing);
71
72#endif
73#endif
74
75          // bug: lm(p) may not divide lm(p * var(i)) in a SCA!
76          if( p_LmDivisibleBy(p, q, currRing) )
77            q = nc_ReduceSpoly(p, q, currRing);
78
79
80#ifdef PDEBUG
81          p_Test(q, currRing);
82#if 0
83          Print("reductum q/p: ");
84          p_Write(q, currRing);
85
86          // Print("With J!\n");
87#endif
88#endif
89
90//          if( q != NULL)
91          q = kNF(J, currQuotient, q, 0, KSTD_NF_NONORM); // in currRing!!!
92
93#ifdef PDEBUG
94          p_Test(q, currRing);
95#if 0
96          Print("NF(J/currQuotient)=> q: "); // !
97          p_Write(q, currRing);
98#endif
99#endif
100          if (q!=NULL)
101          {
102            if (p_IsConstant(q, currRing)) // => return (1)!
103            {
104              p_Delete(&q, currRing);
105              id_Delete(&J, currRing);
106
107              if (K != NULL)
108                id_Delete(&K, currRing);
109
110              ideal Q = idInit(1,1); // ring independent!
111              Q->m[0] = p_One(currRing);
112
113              return(Q);
114            }
115
116//            flag = false;
117
118            // K += q:
119
120            ideal Q = idInit(1,1); // ring independent
121            Q->m[0]=q;
122
123            if( K == NULL )
124              K = Q;
125            else
126            {
127              ideal id_tmp = idSimpleAdd(K, Q); // in currRing
128              id_Delete(&K, currRing);
129              id_Delete(&Q, currRing);
130              K = id_tmp; // K += Q
131            }
132          }
133
134
135        } // if q != NULL
136      } // for all variables
137
138    }
139
140    if (K == NULL) // nothing new: i.e. all elements are two-sided
141      return(J);
142    // now we update GrBasis J with K
143    //    iSize=IDELEMS(J);
144#ifdef PDEBUG
145    idTest(J); // in currRing!
146#if 0
147    Print("J:");
148    idPrint(J);
149    PrintLn();
150#endif // debug
151#endif
152
153
154
155#ifdef PDEBUG
156    idTest(K); // in currRing!
157#if 0
158    Print("+K:");
159    idPrint(K);
160    PrintLn();
161#endif // debug
162#endif
163
164
165    int iSize = idElem(J); // ring independent
166
167    // J += K:
168    ideal id_tmp = idSimpleAdd(J,K); // in currRing
169    id_Delete(&K, currRing); id_Delete(&J, currRing);
170
171#if 1
172    BITSET save_test=test;
173    test|=Sy_bit(OPT_SB_1); // ring independent
174    J = kStd(id_tmp, currQuotient, testHomog, NULL, NULL, 0, iSize); // J = J + K, J - std // in currRing!
175    test = save_test;
176#else
177    J=kStd(id_tmp, currQuotient,testHomog,NULL,NULL,0,0,NULL);
178#endif
179
180    id_Delete(&id_tmp, currRing);
181    idSkipZeroes(J); // ring independent
182
183#ifdef PDEBUG
184    idTest(J); // in currRing!
185#if 0
186    Print("J:");
187    idPrint(J);
188    PrintLn();
189#endif // debug
190#endif
191  } // loop
192}
193
194
195
196
197static ideal idPrepareStd(ideal T, ideal s,  int k)
198{
199  // T is a left SB, without zeros, s is a list with zeros
200#ifdef PDEBUG
201  if (IDELEMS(s)!=IDELEMS(T))
202  {
203    Print("ideals of diff. size!!!");
204  }
205#endif
206  ideal t = idCopy(T);
207  int j,rs=id_RankFreeModule(s, currRing);
208  poly p,q;
209
210  ideal res = idInit(2*idElem(t),1+idElem(t));
211  if (rs == 0)
212  {
213    for (j=0; j<IDELEMS(t); j++)
214    {
215      if (s->m[j]!=NULL) pSetCompP(s->m[j],1);
216      if (t->m[j]!=NULL) pSetCompP(t->m[j],1);
217    }
218    k = si_max(k,1);
219  }
220  for (j=0; j<IDELEMS(t); j++)
221  {
222    if (s->m[j]!=NULL)
223    {
224      p = s->m[j];
225      q = pOne();
226      pSetComp(q,k+1+j);
227      pSetmComp(q);
228#if 0
229      while (pNext(p)) pIter(p);
230      pNext(p) = q;
231#else
232      p = pAdd(p,q);
233      s->m[j] = p;
234#ifdef PDEBUG
235    pTest(p);
236#endif
237#endif
238    }
239  }
240  res = idSimpleAdd(t,s);
241  idDelete(&t);
242  res->rank = 1+idElem(T);
243  return(res);
244}
245
246
247ideal Approx_Step(ideal L)
248{
249  int N=currRing->N;
250  int i,j; // k=syzcomp
251  int flag, flagcnt=0, syzcnt=0;
252  int syzcomp = 0;
253  int k=1; // for ideals not modules
254  ideal I = kStd(L, currQuotient,testHomog,NULL,NULL,0,0,NULL);
255  idSkipZeroes(I);
256  ideal s_I;
257  int idI = idElem(I);
258  ideal trickyQuotient;
259  if (currQuotient !=NULL)
260  {
261    trickyQuotient = idSimpleAdd(currQuotient,I);
262  }
263  else
264    trickyQuotient = I;
265  idSkipZeroes(trickyQuotient);
266  poly *var = (poly *)omAlloc0((N+1)*sizeof(poly));
267  //  poly *W = (poly *)omAlloc0((2*N+1)*sizeof(poly));
268  resolvente S = (resolvente)omAlloc0((N+1)*sizeof(ideal));
269  ideal SI, res;
270  matrix MI;
271  poly x=pOne();
272  var[0]=x;
273  ideal   h2, h3, s_h2, s_h3;
274  poly    p,q,qq;
275  // init vars
276  for (i=1; i<=N; i++ )
277  {
278    x = pOne();
279    pSetExp(x,i,1);
280    pSetm(x);
281    var[i]=pCopy(x);
282  }
283  // init NF's
284  for (i=1; i<=N; i++ )
285  {
286    h2 = idInit(idI,1);
287    flag = 0;
288    for (j=0; j< idI; j++ )
289    {
290      q = pp_Mult_mm(I->m[j],var[i],currRing);
291      q = kNF(I,currQuotient,q,0,0);
292      if (q!=0)
293      {
294    h2->m[j]=pCopy(q);
295    //  p_Shift(&(h2->m[flag]),1, currRing);
296    flag++;
297    pDelete(&q);
298      }
299      else
300    h2->m[j]=0;
301    }
302    // W[1..idElems(I)]
303    if (flag >0)
304    {
305      // compute syzygies with values in I
306      //      idSkipZeroes(h2);
307      //      h2 = idSimpleAdd(h2,I);
308      //      h2->rank=flag+idI+1;
309      idTest(h2);
310      //idShow(h2);
311      ring orig_ring = currRing;
312      ring syz_ring = rAssure_SyzComp(orig_ring, TRUE);
313      syzcomp = 1;
314      rSetSyzComp(syzcomp, syz_ring);
315      if (orig_ring != syz_ring)
316      {
317        rChangeCurrRing(syz_ring);
318        s_h2=idrCopyR_NoSort(h2,orig_ring, syz_ring);
319        //  s_trickyQuotient=idrCopyR_NoSort(trickyQuotient,orig_ring);
320        //  rDebugPrint(syz_ring);
321        s_I=idrCopyR_NoSort(I,orig_ring, syz_ring);
322      }
323      else
324      {
325        s_h2 = h2;
326        s_I  = I;
327        //  s_trickyQuotient=trickyQuotient;
328      }
329      idTest(s_h2);
330      //      idTest(s_trickyQuotient);
331      Print(".proceeding with the variable %d\n",i);
332      s_h3 = idPrepareStd(s_I, s_h2, 1);
333      BITSET save_test=test;
334      test|=Sy_bit(OPT_SB_1);
335      idTest(s_h3);
336      idDelete(&s_h2);
337      s_h2=idCopy(s_h3);
338      idDelete(&s_h3);
339      Print("...computing Syz");
340      s_h3 = kStd(s_h2, currQuotient,(tHomog)FALSE,NULL,NULL,syzcomp,idI);
341      test=save_test;
342      //idShow(s_h3);
343      if (orig_ring != syz_ring)
344      {
345        idDelete(&s_h2);
346        for (j=0; j<IDELEMS(s_h3); j++)
347        {
348          if (s_h3->m[j] != NULL)
349          {
350            if (p_MinComp(s_h3->m[j],syz_ring) > syzcomp) // i.e. it is a syzygy
351              p_Shift(&s_h3->m[j], -syzcomp, currRing);
352            else
353              pDelete(&s_h3->m[j]);
354          }
355        }
356        idSkipZeroes(s_h3);
357        s_h3->rank -= syzcomp;
358        rChangeCurrRing(orig_ring);
359        //  s_h3 = idrMoveR_NoSort(s_h3, syz_ring, orig_ring);
360        s_h3 = idrMoveR_NoSort(s_h3, syz_ring, orig_ring);
361        rDelete(syz_ring);
362      }
363      idTest(s_h3);
364      S[syzcnt]=kStd(s_h3,currQuotient,(tHomog)FALSE,NULL,NULL);
365      syzcnt++;
366      idDelete(&s_h3);
367    } // end if flag >0
368    else
369    {
370      flagcnt++;
371    }
372  }
373  if (flagcnt == N)
374  {
375    Print("the input is a two--sided ideal");
376    return(I);
377  }
378  if (syzcnt >0)
379  {
380    Print("..computing Intersect of %d modules\n",syzcnt);
381    if (syzcnt == 1)
382      SI = S[0];
383    else
384      SI = idMultSect(S, syzcnt);
385    //idShow(SI);
386    MI = id_Module2Matrix(SI,currRing);
387    res= idInit(MATCOLS(MI),1);
388    for (i=1; i<= MATCOLS(MI); i++)
389    {
390      p = NULL;
391      for (j=0; j< idElem(I); j++)
392      {
393        q = pCopy(MATELEM(MI,j+1,i));
394        if (q!=NULL)
395        {
396          q = pMult(q,pCopy(I->m[j]));
397          p = pAdd(p,q);
398        }
399      }
400      res->m[i-1]=p;
401    }
402    Print("final std");
403    res = kStd(res, currQuotient,testHomog,NULL,NULL,0,0,NULL);
404    idSkipZeroes(res);
405    return(res);
406  }
407  else
408  {
409    Print("No syzygies");
410    return(I);
411  }
412}
413
Note: See TracBrowser for help on using the repository browser.