source: git/kernel/ringgb.cc @ 76cfef

spielwiese
Last change on this file since 76cfef was 76cfef, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: fixed #includes in the kernel/ sources ADD: dummy headers in some strange cases
  • Property mode set to 100644
File size: 6.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: ringgb interface
7*/
8//#define HAVE_TAIL_RING
9#define NO_BUCKETS
10
11#include <kernel/mod2.h>
12#include <kernel/kutil.h>
13#include <kernel/structs.h>
14#include <omalloc/omalloc.h>
15#include <polys/polys.h>
16#include <polys/monomials/p_polys.h>
17#include <kernel/ideals.h>
18#include <kernel/febase.h>
19#include <kernel/kstd1.h>
20#include <kernel/khstd.h>
21#include <polys/kbuckets.h>
22#include <polys/weight.h>
23#include <misc/intvec.h>
24#include <kernel/pInline1.h>
25#ifdef HAVE_PLURAL
26#include <polys/nc/nc.h>
27#endif
28
29#include <kernel/ringgb.h>
30
31#ifdef HAVE_RINGS
32poly reduce_poly_fct(poly p, ring r)
33{
34   return kFindZeroPoly(p, r, r);
35}
36
37/*
38 * Returns maximal k, such that
39 * 2^k | n
40 */
41int indexOf2(number n)
42{
43  long test = (long) n;
44  int i = 0;
45  while (test%2 == 0)
46  {
47    i++;
48    test = test / 2;
49  }
50  return i;
51}
52
53/***************************************************************
54 *
55 * Lcm business
56 *
57 ***************************************************************/
58// get m1 = LCM(LM(p1), LM(p2))/LM(p1)
59//     m2 = LCM(LM(p1), LM(p2))/LM(p2)
60BOOLEAN ring2toM_GetLeadTerms(const poly p1, const poly p2, const ring p_r,
61                               poly &m1, poly &m2, const ring m_r)
62{
63  int i;
64  int x;
65  m1 = p_Init(m_r);
66  m2 = p_Init(m_r);
67
68  for (i = p_r->N; i; i--)
69  {
70    x = p_GetExpDiff(p1, p2, i, p_r);
71    if (x > 0)
72    {
73      p_SetExp(m2,i,x, m_r);
74      p_SetExp(m1,i,0, m_r);
75    }
76    else
77    {
78      p_SetExp(m1,i,-x, m_r);
79      p_SetExp(m2,i,0, m_r);
80    }
81  }
82  p_Setm(m1, m_r);
83  p_Setm(m2, m_r);
84  long cp1 = (long) pGetCoeff(p1);
85  long cp2 = (long) pGetCoeff(p2);
86  if (cp1 != 0 && cp2 != 0)
87  {
88    while (cp1%2 == 0 && cp2%2 == 0)
89    {
90      cp1 = cp1 / 2;
91      cp2 = cp2 / 2;
92    }
93  }
94  p_SetCoeff(m1, (number) cp2, m_r);
95  p_SetCoeff(m2, (number) cp1, m_r);
96  return TRUE;
97}
98
99void printPolyMsg(const char * start, poly f, const char * end)
100{
101  PrintS(start);
102  wrp(f);
103  PrintS(end);
104}
105
106poly spolyRing2toM(poly f, poly g, ring r)
107{
108  poly m1 = NULL;
109  poly m2 = NULL;
110  ring2toM_GetLeadTerms(f, g, r, m1, m2, r);
111  // printPolyMsg("spoly: m1=", m1, " | ");
112  // printPolyMsg("m2=", m2, "");
113  // PrintLn();
114  poly sp = pSub(p_Mult_mm(f, m1, r), pp_Mult_mm(g, m2, r));
115  pDelete(&m1);
116  pDelete(&m2);
117  return(sp);
118}
119
120poly ringRedNF (poly f, ideal G, ring r)
121{
122  // If f = 0, then normal form is also 0
123  if (f == NULL) { return NULL; }
124  poly h = NULL;
125  poly g = pCopy(f);
126  int c = 0;
127  while (g != NULL)
128  {
129    Print("%d-step RedNF - g=", c);
130    wrp(g);
131    PrintS(" | h=");
132    wrp(h);
133    PrintLn();
134    g = ringNF(g, G, r);
135    if (g != NULL) {
136      h = pAdd(h, pHead(g));
137      pLmDelete(&g);
138    }
139    c++;
140  }
141  return h;
142}
143
144#endif
145
146#ifdef HAVE_RINGS
147
148/*
149 * Find an index i from G, such that
150 * LT(rside) = x * LT(G[i]) has a solution
151 * or -1 if rside is not in the
152 * ideal of the leading coefficients
153 * of the suitable g from G.
154 */
155int findRingSolver(poly rside, ideal G, ring r)
156{
157  if (rside == NULL) return -1;
158  int i;
159//  int iO2rside = indexOf2(pGetCoeff(rside));
160  for (i = 0; i < IDELEMS(G); i++)
161  {
162    if // (indexOf2(pGetCoeff(G->m[i])) <= iO2rside &&    / should not be necessary any more
163       (p_LmDivisibleBy(G->m[i], rside, r))
164    {
165      return i;
166    }
167  }
168  return -1;
169}
170
171poly plain_spoly(poly f, poly g)
172{
173  number cf = nCopy(pGetCoeff(f)), cg = nCopy(pGetCoeff(g));
174  (void)ksCheckCoeff(&cf, &cg); // gcd and zero divisors
175  poly fm, gm;
176  k_GetLeadTerms(f, g, currRing, fm, gm, currRing);
177  pSetCoeff0(fm, cg);
178  pSetCoeff0(gm, cf);  // and now, m1 * LT(p1) == m2 * LT(p2)
179  poly sp = pSub(ppMult_mm(f, fm), ppMult_mm(g, gm));
180  pDelete(&fm);
181  pDelete(&gm);
182  return(sp);
183}
184
185/*2
186* Generates spoly(0, h) if applicable. Assumes ring in Z/2^n.
187*/
188poly plain_zero_spoly(poly h)
189{
190  poly p = NULL;
191  number gcd = nGcd((number) 0, pGetCoeff(h), currRing);
192  if ((NATNUMBER) gcd > 1)
193  {
194    p = p_Copy(h->next, currRing);
195    p = p_Mult_nn(p, nIntDiv(0, gcd), currRing);
196  }
197  return p;
198}
199
200poly ringNF(poly f, ideal G, ring r)
201{
202  // If f = 0, then normal form is also 0
203  if (f == NULL) { return NULL; }
204  poly tmp = NULL;
205  poly h = pCopy(f);
206  int i = findRingSolver(h, G, r);
207  int c = 1;
208  while (h != NULL && i >= 0) {
209//    Print("%d-step NF - h:", c);
210//    wrp(h);
211//    PrintS(" ");
212//    PrintS("G->m[i]:");
213//    wrp(G->m[i]);
214//    PrintLn();
215    tmp = h;
216    h = plain_spoly(h, G->m[i]);
217    pDelete(&tmp);
218//    PrintS("=> h=");
219//    wrp(h);
220//    PrintLn();
221    i = findRingSolver(h, G, r);
222    c++;
223  }
224  return h;
225}
226
227int testGB(ideal I, ideal GI) {
228  poly f, g, h, nf;
229  int i = 0;
230  int j = 0;
231  PrintS("I included?");
232  for (i = 0; i < IDELEMS(I); i++) {
233    if (ringNF(I->m[i], GI, currRing) != NULL) {
234      PrintS("Not reduced to zero from I: ");
235      wrp(I->m[i]);
236      PrintS(" --> ");
237      wrp(ringNF(I->m[i], GI, currRing));
238      PrintLn();
239      return(0);
240    }
241    PrintS("-");
242  }
243  PrintS(" Yes!\nspoly --> 0?");
244  for (i = 0; i < IDELEMS(GI); i++)
245  {
246    for (j = i + 1; j < IDELEMS(GI); j++)
247    {
248      f = pCopy(GI->m[i]);
249      g = pCopy(GI->m[j]);
250      h = plain_spoly(f, g);
251      nf = ringNF(h, GI, currRing);
252      if (nf != NULL)
253      {
254        PrintS("spoly(");
255        wrp(GI->m[i]);
256        PrintS(", ");
257        wrp(GI->m[j]);
258        PrintS(") = ");
259        wrp(h);
260        PrintS(" --> ");
261        wrp(nf);
262        PrintLn();
263        return(0);
264      }
265      pDelete(&f);
266      pDelete(&g);
267      pDelete(&h);
268      pDelete(&nf);
269      PrintS("-");
270    }
271  }
272  if (!(rField_is_Domain()))
273  {
274    PrintS(" Yes!\nzero-spoly --> 0?");
275    for (i = 0; i < IDELEMS(GI); i++)
276    {
277      f = plain_zero_spoly(GI->m[i]);
278      nf = ringNF(f, GI, currRing);
279      if (nf != NULL) {
280        PrintS("spoly(");
281        wrp(GI->m[i]);
282        PrintS(", ");
283        wrp(0);
284        PrintS(") = ");
285        wrp(h);
286        PrintS(" --> ");
287        wrp(nf);
288        PrintLn();
289        return(0);
290      }
291      pDelete(&f);
292      pDelete(&nf);
293      PrintS("-");
294    }
295  }
296  PrintS(" Yes!");
297  PrintLn();
298  return(1);
299}
300
301#endif
Note: See TracBrowser for help on using the repository browser.