source: git/kernel/ringgb.cc @ fbc7cb

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