source: git/kernel/nc.cc @ fbc7cb

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