shiftop.cc
Go to the documentation of this file.
1 #include "shiftop.h"
2 
3 #ifdef HAVE_SHIFTBBA
4 
5 #include "templates/p_MemCopy.h"
6 #include "monomials/p_polys.h"
7 #include "polys/simpleideals.h"
8 
9 /* #define SHIFT_MULT_DEBUG */
10 
11 /* enable compat mode until the user interface is updated to support xy instead of x(1)*y(2)
12  * NOTE: it already works, but all tests and the libraries need to be updated first
13  * -> wait until the new interface is released
14 */
15 #define SHIFT_MULT_COMPAT_MODE
16 
17 #ifdef SHIFT_MULT_DEBUG
18 #include "../kernel/polys.h"
19 #endif
20 
21 poly shift_pp_Mult_mm(poly p, const poly m, const ring ri)
22 {
23 #ifdef SHIFT_MULT_DEBUG
24  PrintLn(); PrintS("shift_pp_Mult_mm: ("); p_wrp(p, ri, ri); PrintS(") * "); p_wrp(m, ri, ri);
25 #endif
26 
27  p_Test(p, ri);
28  p_LmTest(m, ri);
29  if (p == NULL)
30  {
31  return NULL;
32  }
33 
34  int lV = ri->isLPring;
35  poly _m = m; // temp hack because m is const
36 #ifdef SHIFT_MULT_COMPAT_MODE
37  _m = p_Copy(_m, ri);
38  p_mLPunshift(_m, ri);
39  p = p_Copy(p, ri);
40  poly pCopyHead = p; // used to delete p later
41  p_LPunshift(p, ri);
42 #else
43  assume(p_mFirstVblock(_m, ri) <= 1);
44  assume(p_FirstVblock(p, ri) <= 1); // TODO check that each block is <=1
45 #endif
46  // at this point _m and p are shifted to 1
47 
48  spolyrec rp;
49  poly q = &rp; // we use p for iterating and q for the result
50  number mCoeff = pGetCoeff(_m);
51  omBin bin = ri->PolyBin;
52  pAssume(!n_IsZero(mCoeff, ri->cf));
53  pAssume1(p_GetComp(m, ri) == 0 || p_MaxComp(p, ri) == 0);
54 
55  int *mExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
56  p_GetExpV(_m,mExpV,ri);
57  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
58  int *pExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
59  do
60  {
61  p_AllocBin(pNext(q), bin, ri);
62  pIter(q);
63  pSetCoeff0(q, n_Mult(mCoeff, pGetCoeff(p), ri->cf));
64 
65  p_GetExpV(p, pExpV, ri);
66  p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
67  p_MemCopy_LengthGeneral(q->exp, p->exp, ri->ExpL_Size); // otherwise q is not initialized correctly
68  p_SetExpV(q, pExpV, ri);
69 
70  pIter(p);
71  }
72  while (p != NULL);
73  omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
74  omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
75  pNext(q) = NULL;
76 #ifdef SHIFT_MULT_COMPAT_MODE
77  p_Delete(&_m, ri); // in this case we copied _m before
78  p_Delete(&pCopyHead, ri); // in this case we copied p before
79 #endif
80 #ifdef SHIFT_MULT_DEBUG
81  PrintLn(); PrintS("shift_pp_Mult_mm result: "); p_wrp(pNext(&rp), ri, ri); PrintLn();
82 #endif
83  p_Test(pNext(&rp), ri);
84  return pNext(&rp);
85 }
86 
87 // destroys p
88 poly shift_p_Mult_mm(poly p, const poly m, const ring ri)
89 {
90 #ifdef SHIFT_MULT_DEBUG
91  PrintLn(); PrintS("shift_p_Mult_mm: ("); p_wrp(p, ri, ri); PrintS(") * "); p_wrp(m, ri, ri);
92 #endif
93 
94  p_Test(p, ri);
95  p_LmTest(m, ri);
96  pAssume(m != NULL);
97  assume(p!=NULL);
98 
99  int lV = ri->isLPring;
100  poly _m = m; // temp hack because m is const
101 #ifdef SHIFT_MULT_COMPAT_MODE
102  _m = p_Copy(_m, ri);
103  p_mLPunshift(_m, ri);
104  p_LPunshift(p, ri);
105 #else
106  assume(p_mFirstVblock(_m, ri) <= 1);
107  assume(p_FirstVblock(p, ri) <= 1); // TODO check that each block is <=1
108 #endif
109  // at this point _m and p are shifted to 1
110 
111  poly q = p; // we use p for iterating and q for the result
112  number mCoeff = pGetCoeff(_m);
113  number pCoeff;
114  pAssume(!n_IsZero(mCoeff, ri->cf));
115 
116  int *mExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
117  p_GetExpV(_m,mExpV,ri);
118  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
119  int *pExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
120  while (p != NULL)
121  {
122  pCoeff = pGetCoeff(p);
123  pSetCoeff0(p, n_Mult(mCoeff, pCoeff, ri->cf));
124  n_Delete(&pCoeff, ri->cf); // delete the old coeff
125 
126  p_GetExpV(p,pExpV,ri);
127  p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
128  p_SetExpV(p, pExpV, ri);
129 
130  pIter(p);
131  }
132  omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
133  omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
134 #ifdef SHIFT_MULT_COMPAT_MODE
135  p_Delete(&_m, ri); // in this case we copied _m before
136 #endif
137 #ifdef SHIFT_MULT_DEBUG
138  PrintLn(); PrintS("shift_p_Mult_mm result: "); p_wrp(q, ri, ri); PrintLn();
139 #endif
140  p_Test(q, ri);
141  return q;
142 }
143 
144 poly shift_pp_mm_Mult(poly p, const poly m, const ring ri)
145 {
146 #ifdef SHIFT_MULT_DEBUG
147  PrintLn(); PrintS("shift_pp_mm_Mult: "); p_wrp(m, ri, ri); PrintS(" * ("); p_wrp(p, ri, ri); PrintS(")");
148 #endif
149 
150  p_Test(p, ri);
151  p_LmTest(m, ri);
152  if (p == NULL)
153  {
154  return NULL;
155  }
156 
157  int lV = ri->isLPring;
158  poly _m = m; // temp hack because m is const
159 #ifdef SHIFT_MULT_COMPAT_MODE
160  _m = p_Copy(_m, ri);
161  p_mLPunshift(_m, ri);
162  p = p_Copy(p, ri);
163  poly pCopyHead = p; // used to delete p later
164  p_LPunshift(p, ri);
165 #else
166  assume(p_mFirstVblock(_m, ri) <= 1);
167  assume(p_FirstVblock(p, ri) <= 1); // TODO check that each block is <=1
168 #endif
169  // at this point _m and p are shifted to 1
170 
171  spolyrec rp;
172  poly q = &rp; // we use p for iterating and q for the result
173  number mCoeff = pGetCoeff(_m);
174  omBin bin = ri->PolyBin;
175  pAssume(!n_IsZero(mCoeff, ri->cf));
176  pAssume1(p_GetComp(m, ri) == 0 || p_MaxComp(p, ri) == 0);
177 
178  int *mExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
179  p_GetExpV(_m,mExpV,ri);
180  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
181  int *pExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
182  do
183  {
184  p_AllocBin(pNext(q), bin, ri);
185  pIter(q);
186  pSetCoeff0(q, n_Mult(mCoeff, pGetCoeff(p), ri->cf));
187 
188  p_GetExpV(p, pExpV, ri);
189  p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
190  p_MemCopy_LengthGeneral(q->exp, p->exp, ri->ExpL_Size); // otherwise q is not initialized correctly
191  p_SetExpV(q, pExpV, ri);
192 
193  pIter(p);
194  }
195  while (p != NULL);
196  omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
197  omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
198  pNext(q) = NULL;
199 #ifdef SHIFT_MULT_COMPAT_MODE
200  p_Delete(&_m, ri); // in this case we copied _m before
201  p_Delete(&pCopyHead, ri); // in this case we copied p before
202 #endif
203 #ifdef SHIFT_MULT_DEBUG
204  PrintLn(); PrintS("shift_pp_mm_Mult result: "); p_wrp(pNext(&rp), ri, ri); PrintLn();
205 #endif
206  p_Test(pNext(&rp), ri);
207  return pNext(&rp);
208 }
209 
210 // destroys p
211 poly shift_p_mm_Mult(poly p, const poly m, const ring ri)
212 {
213 #ifdef SHIFT_MULT_DEBUG
214  PrintLn(); PrintS("shift_p_mm_Mult: "); p_wrp(m, ri, ri); PrintS(" * ("); p_wrp(p, ri, ri); PrintS(")");
215 #endif
216 
217  p_Test(p, ri);
218  p_LmTest(m, ri);
219  pAssume(m != NULL);
220  assume(p!=NULL);
221 
222  int lV = ri->isLPring;
223  poly _m = m; // temp hack because m is const
224 #ifdef SHIFT_MULT_COMPAT_MODE
225  _m = p_Copy(_m, ri);
226  p_mLPunshift(_m, ri);
227  p_LPunshift(p, ri);
228 #else
229  assume(p_mFirstVblock(_m, ri) <= 1);
230  assume(p_FirstVblock(p, ri) <= 1); // TODO check that each block is <=1
231 #endif
232  // at this point _m and p are shifted to 1
233 
234  poly q = p; // we use p for iterating and q for the result
235  number mCoeff = pGetCoeff(_m);
236  number pCoeff;
237  pAssume(!n_IsZero(mCoeff, ri->cf));
238 
239  int *mExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
240  p_GetExpV(_m,mExpV,ri);
241  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
242  int *pExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
243  while (p != NULL)
244  {
245  pCoeff = pGetCoeff(p);
246  pSetCoeff0(p, n_Mult(mCoeff, pCoeff, ri->cf));
247  n_Delete(&pCoeff, ri->cf); // delete the old coeff
248 
249  p_GetExpV(p,pExpV,ri);
250  p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
251  p_SetExpV(p, pExpV, ri);
252 
253  pIter(p);
254  }
255  omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
256  omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
257 #ifdef SHIFT_MULT_COMPAT_MODE
258  p_Delete(&_m, ri); // in this case we copied _m before
259 #endif
260 #ifdef SHIFT_MULT_DEBUG
261  PrintLn(); PrintS("shift_p_mm_Mult result: "); p_wrp(q, ri, ri); PrintLn();
262 #endif
263  p_Test(q, ri);
264  return q;
265 }
266 
267 // p - m*q destroys p
268 poly shift_p_Minus_mm_Mult_qq(poly p, poly m, poly q, int& Shorter, const poly spNoether, const ring ri) {
269 #ifdef SHIFT_MULT_DEBUG
270  PrintLn(); PrintS("shift_p_Minus_mm_Mult_qq: "); p_wrp(p, ri, ri); PrintS(" - "); p_wrp(m, ri, ri); PrintS(" * "); p_wrp(q, ri, ri);
271 #endif
272 
273  Shorter = pLength(p) + pLength(q);
274 
275  poly qq = p_Add_q(p, shift_pp_mm_Mult(q, p_Neg(p_Copy(m, ri), ri), ri), ri);
276 
277 #ifdef SHIFT_MULT_DEBUG
278  PrintLn(); PrintS("shift_p_Minus_mm_Mult_qq result: "); p_wrp(qq, ri, ri); PrintLn();
279 #endif
280  Shorter -= pLength(qq);
281  return qq;
282 }
283 
284 // Unsupported Operation STUBs
285 poly shift_pp_Mult_mm_Noether_STUB(poly p, const poly m, const poly spNoether, int &ll, const ring ri) {
286  PrintLn(); WarnS("pp_Mult_mm_Noether is not supported yet by Letterplace. Ignoring spNoether and using pp_Mult_mm. This might lead to unexpected behavior.");
287 
288  int pLen = 0;
289  if (ll >= 0)
290  {
291  pLen = pLength(p);
292  }
293 
294  p = shift_pp_Mult_mm(p, m, ri);
295 
296  if (ll >= 0)
297  {
298  ll = pLen - pLength(p);
299  }
300  else
301  {
302  ll = pLength(p);
303  }
304 
305  return p;
306 }
307 
308 
309 poly shift_pp_Mult_Coeff_mm_DivSelectMult_STUB(poly p,const poly m, const poly a, const poly b, int &shorter,const ring r) {
310  PrintLn(); WarnS("pp_Mult_Coeff_mm_DivSelectMult is not supported yet by Letterplace. This might lead to unexpected behavior.");
311  return NULL;
312 }
313 
314 poly shift_pp_Mult_Coeff_mm_DivSelect_STUB(poly p, const poly m, int &shorter, const ring r) {
315  PrintLn(); WarnS("pp_Mult_Coeff_mm_DivSelect is not supported yet by Letterplace. This might lead to unexpected behavior.");
316  return NULL;
317 }
318 
319 // auxiliary
320 
321 // unshifts the monomial m
322 void p_mLPunshift(poly m, const ring ri)
323 {
324  if (m == NULL || p_LmIsConstantComp(m,ri)) return;
325 
326  int lV = ri->isLPring;
327 
328  int shift = p_mFirstVblock(m, ri) - 1;
329 
330  if (shift == 0) return;
331 
332  int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
333  int *s=(int *)omAlloc0((ri->N+1)*sizeof(int));
334  p_GetExpV(m, e, ri);
335 
336  int expVoffset = shift*lV;
337  for (int i = 1 + expVoffset; i <= ri->N; i++)
338  {
339  assume(e[i] <= 1);
340  s[i - expVoffset] = e[i];
341  }
342  p_SetExpV(m,s,ri);
343  omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
344  omFreeSize((ADDRESS) s, (ri->N+1)*sizeof(int));
345 }
346 
347 // unshifts the polynomial p, note: the ordering can be destroyed if the shifts for the monomials are not equal
348 void p_LPunshift(poly p, const ring ri)
349 {
350  while (p!=NULL)
351  {
352  p_mLPunshift(p, ri);
353  pIter(p);
354  }
355 }
356 
357 void p_mLPshift(poly m, int sh, const ring ri)
358 {
359  if (sh == 0 || m == NULL || p_LmIsConstantComp(m,ri)) return;
360 
361  int lV = ri->isLPring;
362 
363  assume(p_mFirstVblock(m,ri) + sh >= 1);
364  assume(p_mLastVblock(m,ri) + sh <= ri->N/lV);
365 
366  int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
367  int *s=(int *)omAlloc0((ri->N+1)*sizeof(int));
368  p_GetExpV(m,e,ri);
369 
370  for (int i = ri->N - sh*lV; i > 0; i--)
371  {
372  assume(e[i]<=1);
373  if (e[i]==1)
374  {
375  s[i + (sh*lV)] = e[i]; /* actually 1 */
376  }
377  }
378  p_SetExpV(m,s,ri);
379  omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
380  omFreeSize((ADDRESS) s, (ri->N+1)*sizeof(int));
381 }
382 
383 void p_LPshift(poly p, int sh, const ring ri)
384 {
385  if (sh == 0) return;
386 
387  while (p!=NULL)
388  {
389  p_mLPshift(p, sh, ri);
390  pIter(p);
391  }
392 }
393 
394 /* returns the number of maximal block */
395 /* appearing among the monomials of p */
396 /* the 0th block is the 1st one */
397 int p_LastVblock(poly p, const ring r)
398 {
399  poly q = p;
400  int ans = 0;
401  while (q!=NULL)
402  {
403  int ansnew = p_mLastVblock(q, r);
404  ans = si_max(ans,ansnew);
405  pIter(q);
406  }
407  return(ans);
408 }
409 
410 /* for a monomial p, returns the number of the last block */
411 /* where a nonzero exponent is sitting */
412 int p_mLastVblock(poly p, const ring ri)
413 {
414  if (p == NULL || p_LmIsConstantComp(p,ri))
415  {
416  return(0);
417  }
418 
419  int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
420  p_GetExpV(p,e,ri);
421  int b = p_mLastVblock(p, e, ri);
422  omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
423  return b;
424 }
425 
426 /* for a monomial p with exponent vector expV, returns the number of the last block */
427 /* where a nonzero exponent is sitting */
428 int p_mLastVblock(poly p, int *expV, const ring ri)
429 {
430  if (p == NULL || p_LmIsConstantComp(p,ri))
431  {
432  return(0);
433  }
434 
435  int lV = ri->isLPring;
436  int j,b;
437  j = ri->N;
438  while ( (!expV[j]) && (j>=1) ) j--;
439  assume(j>0);
440  b = (int)((j+lV-1)/lV); /* the number of the block, >=1 */
441  return b;
442 }
443 
444 /* returns the number of maximal block */
445 /* appearing among the monomials of p */
446 /* the 0th block is the 1st one */
447 int p_FirstVblock(poly p, const ring r)
448 {
449  if (p == NULL) {
450  return 0;
451  }
452 
453  poly q = p;
454  int ans = p_mFirstVblock(q, r);
455  while (q!=NULL)
456  {
457  int ansnew = p_mFirstVblock(q, r);
458  if (ansnew > 0) { // don't count constants
459  ans = si_min(ans,ansnew);
460  }
461  pIter(q);
462  }
463  /* do not need to delete q */
464  return(ans);
465 }
466 
467 /* for a monomial p, returns the number of the first block */
468 /* where a nonzero exponent is sitting */
469 int p_mFirstVblock(poly p, const ring ri)
470 {
471  if (p == NULL || p_LmIsConstantComp(p,ri))
472  {
473  return(0);
474  }
475 
476  int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
477  p_GetExpV(p,e,ri);
478  int b = p_mFirstVblock(p, e, ri);
479  omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
480  return b;
481 }
482 
483 /* for a monomial p with exponent vector expV, returns the number of the first block */
484 /* where a nonzero exponent is sitting */
485 int p_mFirstVblock(poly p, int *expV, const ring ri)
486 {
487  if (p == NULL || p_LmIsConstantComp(p,ri))
488  {
489  return(0);
490  }
491 
492  int lV = ri->isLPring;
493  int j,b;
494  j = 1;
495  while ( (!expV[j]) && (j<=ri->N-1) ) j++;
496  assume(j <= ri->N);
497  b = (int)(j+lV-1)/lV; /* the number of the block, 1<= b <= r->N */
498  return b;
499 }
500 
501 // appends m2ExpV to m1ExpV, also adds their components (one of them is always zero)
502 void p_LPExpVappend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri) {
503 #ifdef SHIFT_MULT_DEBUG
504  PrintLn(); PrintS("Append");
505  PrintLn(); WriteLPExpV(m1ExpV, ri);
506  PrintLn(); WriteLPExpV(m2ExpV, ri);
507 #endif
508  if (m1Length + m2Length > ri->N)
509  {
510  WarnS("letterplace degree bound too low for this multiplication");
511  }
512  for (int i = 1 + m1Length; i < 1 + m1Length + m2Length; ++i)
513  {
514  assume(m2ExpV[i - m1Length] <= 1);
515  m1ExpV[i] = m2ExpV[i - m1Length];
516  }
517 
518  assume(m1ExpV[0] == 0 || m2ExpV[0] == 0); // one component should be zero (otherwise this doesn't make any sense)
519  m1ExpV[0] += m2ExpV[0]; // as in the commutative variant (they use MemAdd)
520 #ifdef SHIFT_MULT_DEBUG
521  PrintLn(); WriteLPExpV(m1ExpV, ri);
522 #endif
523 }
524 
525 // prepends m2ExpV to m1ExpV, also adds their components (one of them is always zero)
526 void p_LPExpVprepend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri)
527 {
528 #ifdef SHIFT_MULT_DEBUG
529  PrintLn(); PrintS("Prepend");
530  PrintLn(); WriteLPExpV(m1ExpV, ri);
531  PrintLn(); WriteLPExpV(m2ExpV, ri);
532 #endif
533  if (m1Length + m2Length > ri->N)
534  {
535  WarnS("letterplace degree bound too low for this multiplication");
536  }
537 
538  // shift m1 by m2Length
539  for (int i = m2Length + m1Length; i >= 1 + m2Length; --i)
540  {
541  m1ExpV[i] = m1ExpV[i - m2Length];
542  }
543 
544  // write m2 to m1
545  for (int i = 1; i < 1 + m2Length; ++i)
546  {
547  assume(m2ExpV[i] <= 1);
548  m1ExpV[i] = m2ExpV[i];
549  }
550 
551  assume(m1ExpV[0] == 0 || m2ExpV[0] == 0); // one component should be zero (otherwise this doesn't make any sense)
552  m1ExpV[0] += m2ExpV[0]; // as in the commutative variant (they use MemAdd)
553 #ifdef SHIFT_MULT_DEBUG
554  PrintLn(); WriteLPExpV(m1ExpV, ri);
555 #endif
556 }
557 
558 void WriteLPExpV(int *expV, ring ri)
559 {
560  char *s = LPExpVString(expV, ri);
561  PrintS(s);
562  omFree(s);
563 }
564 
565 char* LPExpVString(int *expV, ring ri)
566 {
567  StringSetS("");
568  for (int i = 0; i <= ri->N; ++i)
569  {
570  StringAppend("%d", expV[i]);
571  if (i == 0)
572  {
573  StringAppendS("| ");
574  }
575  if (i % ri->isLPring == 0 && i != ri->N)
576  {
577  StringAppendS(" ");
578  }
579  }
580  return StringEndS();
581 }
582 
583 // splits a frame (e.g. x(1)*y(5)) m1 into m1 and m2 (e.g. m1=x(1) and m2=y(1))
584 // according to p which is inside the frame
585 void k_SplitFrame(poly &m1, poly &m2, int at, const ring r)
586 {
587  int lV = r->isLPring;
588 
589  number m1Coeff = pGetCoeff(m1);
590 
591  int hole = lV * at;
592  m2 = p_GetExp_k_n(m1, 1, hole, r);
593  m1 = p_GetExp_k_n(m1, hole, r->N, r);
594 
595  p_mLPunshift(m2, r);
596  p_SetCoeff(m1, m1Coeff, r);
597 
598  assume(p_FirstVblock(m1,r) <= 1);
599  assume(p_FirstVblock(m2,r) <= 1);
600 }
601 
602 /* tests whether each polynomial of an ideal I lies in in V */
603 int id_IsInV(ideal I, const ring r)
604 {
605  int i;
606  int s = IDELEMS(I)-1;
607  for(i = 0; i <= s; i++)
608  {
609  if ( !p_IsInV(I->m[i], r) )
610  {
611  return(0);
612  }
613  }
614  return(1);
615 }
616 
617 /* tests whether the whole polynomial p in in V */
618 int p_IsInV(poly p, const ring r)
619 {
620  poly q = p;
621  while (q!=NULL)
622  {
623  if ( !p_mIsInV(q, r) )
624  {
625  return(0);
626  }
627  q = pNext(q);
628  }
629  return(1);
630 }
631 
632 /* there should be two routines: */
633 /* 1. test place-squarefreeness: in homog this suffices: isInV */
634 /* 2. test the presence of a hole -> in the tail??? */
635 
636 int p_mIsInV(poly p, const ring r)
637 {
638  int lV = r->isLPring;
639  /* investigate only the leading monomial of p in currRing */
640  if ( p_Totaldegree(p, r)==0 ) return(1);
641  /* returns 1 iff p is in V */
642  /* that is in each block up to a certain one there is only one nonzero exponent */
643  /* lV = the length of V = the number of orig vars */
644  int *e = (int *)omAlloc((r->N+1)*sizeof(int));
645  int b = (int)((r->N+lV-1)/lV); /* the number of blocks */
646  //int b = (int)(currRing->N)/lV;
647  int *B = (int *)omAlloc0((b+1)*sizeof(int)); /* the num of elements in a block */
648  p_GetExpV(p,e,r);
649  int i,j;
650  for (j=1; j<=b; j++)
651  {
652  /* we go through all the vars */
653  /* by blocks in lV vars */
654  for (i=(j-1)*lV + 1; i<= j*lV; i++)
655  {
656  if (e[i]) B[j] = B[j]+1;
657  }
658  }
659  // j = b;
660  // while ( (!B[j]) && (j>=1)) j--;
661  for (j=b; j>=1; j--)
662  {
663  if (B[j]!=0) break;
664  }
665  /* do not need e anymore */
666  omFreeSize((ADDRESS) e, (r->N+1)*sizeof(int));
667 
668  if (j==0) goto ret_true;
669 // {
670 // /* it is a zero exp vector, which is in V */
671 // freeT(B, b);
672 // return(1);
673 // }
674  /* now B[j] != 0 and we test place-squarefreeness */
675  for (; j>=1; j--)
676  {
677  if (B[j]!=1)
678  {
679  omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
680  return(0);
681  }
682  }
683  ret_true:
684  omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
685  return(1);
686 }
687 
688 BOOLEAN p_LPDivisibleBy(poly a, poly b, const ring r)
689 {
690  pIfThen1(b!=NULL, p_LmCheckPolyRing1(b, r));
691  pIfThen1(a!=NULL, p_LmCheckPolyRing1(a, r));
692 
693  if (b == NULL) return TRUE;
694  if (a != NULL && (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r)))
695  return _p_LPLmDivisibleByNoComp(a,b,r);
696  return FALSE;
697 }
698 
699 BOOLEAN p_LPLmDivisibleBy(poly a, poly b, const ring r)
700 {
701  p_LmCheckPolyRing1(b, r);
702  pIfThen1(a != NULL, p_LmCheckPolyRing1(b, r));
703  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
704  return _p_LPLmDivisibleByNoComp(a, b, r);
705  return FALSE;
706 }
707 
708 BOOLEAN _p_LPLmDivisibleByNoComp(poly a, poly b, const ring r)
709 {
710  if(p_LmIsConstantComp(a, r))
711  return TRUE;
712 #ifdef SHIFT_MULT_COMPAT_MODE
713  a = p_Head(a, r);
714  p_mLPunshift(a, r);
715  b = p_Head(b, r);
716  p_mLPunshift(b, r);
717 #endif
718  int i = (r->N / r->isLPring) - p_LastVblock(a, r);
719  do {
720  int j = r->N - (i * r->isLPring);
721  bool divisible = true;
722  do
723  {
724  if (p_GetExp(a, j, r) > p_GetExp(b, j + (i * r->isLPring), r))
725  {
726  divisible = false;
727  break;
728  }
729  j--;
730  }
731  while (j);
732  if (divisible) return TRUE;
733  i--;
734  }
735  while (i > -1);
736 #ifdef SHIFT_MULT_COMPAT_MODE
737  p_Delete(&a, r);
738  p_Delete(&b, r);
739 #endif
740  return FALSE;
741 }
742 
743 poly p_LPVarAt(poly p, int pos, const ring r)
744 {
745  if (p == NULL || pos < 1 || pos > (r->N / r->isLPring)) return NULL;
746  poly v = p_One(r);
747  for (int i = (pos-1) * r->isLPring + 1; i <= pos * r->isLPring; i++) {
748  if (p_GetExp(p, i, r)) {
749  p_SetExp(v, i - (pos-1) * r->isLPring, 1, r);
750  return v;
751  }
752  }
753  return v;
754 }
755 
756 /// substitute weights from orderings a,wp,Wp
757 /// by d copies of it at position p
758 static BOOLEAN freeAlgebra_weights(const ring old_ring, ring new_ring, int p, int d)
759 {
760  omFree(new_ring->wvhdl[p]);
761  int *w=(int*)omAlloc(new_ring->N*sizeof(int));
762  for(int b=0;b<d;b++)
763  {
764  for(int i=old_ring->N-1;i>=0;i--)
765  {
766  if (old_ring->wvhdl[p][i]<-0) return TRUE;
767  w[b*old_ring->N+i]=old_ring->wvhdl[p][i];
768  }
769  }
770  new_ring->wvhdl[p]=w;
771  new_ring->block1[p]=new_ring->N;
772  return FALSE;
773 }
774 
775 ring freeAlgebra(ring r, int d)
776 {
777  ring R=rCopy0(r);
778  int p;
779  if((r->order[0]==ringorder_C)
780  ||(r->order[0]==ringorder_c))
781  p=1;
782  else
783  p=0;
784  // create R->N
785  R->N=r->N*d;
786  R->isLPring=r->N;
787  // create R->order
788  BOOLEAN has_order_a=FALSE;
789  while (r->order[p]==ringorder_a)
790  {
791  if (freeAlgebra_weights(r,R,p,d))
792  {
793  WerrorS("weights must be positive");
794  return NULL;
795  }
796  has_order_a=TRUE;
797  p++;
798  }
799  R->block1[p]=R->N; /* only dp,Dp,wp,Wp; will be discarded for lp*/
800  switch(r->order[p])
801  {
802  case ringorder_dp:
803  case ringorder_Dp:
804  break;
805  case ringorder_wp:
806  case ringorder_Wp:
807  if (freeAlgebra_weights(r,R,p,d))
808  {
809  WerrorS("weights must be positive");
810  return NULL;
811  }
812  break;
813  case ringorder_lp:
814  case ringorder_rp:
815  {
816  if(has_order_a)
817  {
818  WerrorS("ordering (a(..),lp/rp not implemented for LP-rings");
819  return NULL;
820  }
821  int ** wvhdl=(int**)omAlloc0((r->N+3)*sizeof(int*));
822  rRingOrder_t* ord=(rRingOrder_t*)omAlloc0((r->N+3)*sizeof(rRingOrder_t));
823  int* blk0=(int*)omAlloc0((r->N+3)*sizeof(int));
824  int* blk1=(int*)omAlloc0((r->N+3)*sizeof(int));
825  omFree(R->wvhdl); R->wvhdl=wvhdl;
826  omFree(R->order); R->order=ord;
827  omFree(R->block0); R->block0=blk0;
828  omFree(R->block1); R->block1=blk1;
829  for(int i=0;i<r->N;i++)
830  {
831  ord[i+p]=ringorder_a;
832  //Print("entry:%d->a\n",i+p);
833  blk0[i+p]=1;
834  blk1[i+p]=R->N;
835  wvhdl[i+p]=(int*)omAlloc0(R->N*sizeof(int));
836  for(int j=0;j<d;j++)
837  {
838  assume(j*r->N+i<R->N);
839  if (r->order[p]==ringorder_lp)
840  wvhdl[i+p][j*r->N+i]=1;
841  else
842  wvhdl[i+p][(j+1)*r->N-i-1]=1;
843  }
844  }
845  ord[r->N+p]=r->order[p]; /* lp or rp */
846  //Print("entry:%d->lp\n",r->N+p);
847  blk0[r->N+p]=1;
848  blk1[r->N+p]=R->N;
849  // copy component order
850  if (p==1) ord[0]=r->order[0];
851  else if (p==0) ord[r->N+1]=r->order[1];
852  else
853  { // should never happen:
854  WerrorS("ordering not implemented for LP-rings");
855  return NULL;
856  }
857  //if (p==1) PrintS("entry:0 ->c/C\n");
858  //else if (p==0) Print("entry:%d ->c/C\n",r->N+1);
859  break;
860  }
861  default: WerrorS("ordering not implemented for LP-rings");
862  return NULL;
863  }
864  // create R->names
865  char **names=(char**)omAlloc(R->N*sizeof(char*));
866  for(int b=0;b<d;b++)
867  {
868  for(int i=r->N-1;i>=0;i--)
869  names[b*r->N+i]=omStrDup(r->names[i]);
870  }
871  for(int i=r->N-1;i>=0;i--) omFree(R->names[i]);
872  omFree(R->names);
873  R->names=names;
874 
875  rComplete(R,TRUE);
876  return R;
877 }
878 #endif
const CanonicalForm int s
Definition: facAbsFact.cc:55
int j
Definition: facHensel.cc:105
int p_mLastVblock(poly p, const ring ri)
Definition: shiftop.cc:412
static BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
Definition: p_polys.h:961
omBin_t * omBin
Definition: omStructs.h:12
void PrintLn()
Definition: reporter.cc:310
poly shift_pp_Mult_Coeff_mm_DivSelectMult_STUB(poly p, const poly m, const poly a, const poly b, int &shorter, const ring r)
Definition: shiftop.cc:309
static int si_min(const int a, const int b)
Definition: auxiliary.h:139
poly shift_p_mm_Mult(poly p, const poly m, const ring ri)
Definition: shiftop.cc:211
#define FALSE
Definition: auxiliary.h:94
int id_IsInV(ideal I, const ring r)
Definition: shiftop.cc:603
#define pAssume(cond)
Definition: monomials.h:90
#define p_GetComp(p, r)
Definition: monomials.h:64
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1455
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int p_mIsInV(poly p, const ring r)
Definition: shiftop.cc:636
int p_mFirstVblock(poly p, const ring ri)
Definition: shiftop.cc:469
static BOOLEAN freeAlgebra_weights(const ring old_ring, ring new_ring, int p, int d)
substitute weights from orderings a,wp,Wp by d copies of it at position p
Definition: shiftop.cc:758
#define TRUE
Definition: auxiliary.h:98
#define p_AllocBin(p, bin, r)
Definition: monomials.h:248
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1442
void * ADDRESS
Definition: auxiliary.h:133
poly shift_pp_Mult_Coeff_mm_DivSelect_STUB(poly p, const poly m, int &shorter, const ring r)
Definition: shiftop.cc:314
poly shift_pp_mm_Mult(poly p, const poly m, const ring ri)
Definition: shiftop.cc:144
void WerrorS(const char *s)
Definition: feFopen.cc:24
BOOLEAN p_LPLmDivisibleBy(poly a, poly b, const ring r)
Definition: shiftop.cc:699
char * StringEndS()
Definition: reporter.cc:151
void p_LPExpVappend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri)
Definition: shiftop.cc:502
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:44
#define WarnS
Definition: emacs.cc:78
poly shift_pp_Mult_mm_Noether_STUB(poly p, const poly m, const poly spNoether, int &ll, const ring ri)
Definition: shiftop.cc:285
void p_LPExpVprepend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri)
Definition: shiftop.cc:526
#define omAlloc(size)
Definition: omAllocDecl.h:210
char * LPExpVString(int *expV, ring ri)
Definition: shiftop.cc:565
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:411
static void p_SetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1470
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:811
#define pIter(p)
Definition: monomials.h:37
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of &#39;a&#39; and &#39;b&#39;, i.e., a*b
Definition: coeffs.h:636
int p_IsInV(poly p, const ring r)
Definition: shiftop.cc:618
poly shift_p_Mult_mm(poly p, const poly m, const ring ri)
Definition: shiftop.cc:88
CanonicalForm b
Definition: cfModGcd.cc:4044
static poly p_Head(poly p, const ring r)
Definition: p_polys.h:824
void p_LPshift(poly p, int sh, const ring ri)
Definition: shiftop.cc:383
poly shift_p_Minus_mm_Mult_qq(poly p, poly m, poly q, int &Shorter, const poly spNoether, const ring ri)
Definition: shiftop.cc:268
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:48
poly p_One(const ring r)
Definition: p_polys.cc:1303
BOOLEAN rComplete(ring r, int force)
this needs to be called whenever a new ring is created: new fields in ring are created (like VarOffse...
Definition: ring.cc:3394
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
Definition: p_polys.h:468
void WriteLPExpV(int *expV, ring ri)
Definition: shiftop.cc:558
#define omFree(addr)
Definition: omAllocDecl.h:261
#define assume(x)
Definition: mod2.h:390
void StringSetS(const char *st)
Definition: reporter.cc:128
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1365
void StringAppendS(const char *st)
Definition: reporter.cc:107
rRingOrder_t
order stuff
Definition: ring.h:67
int m
Definition: cfEzgcd.cc:121
BOOLEAN _p_LPLmDivisibleByNoComp(poly a, poly b, const ring r)
Definition: shiftop.cc:708
static int si_max(const int a, const int b)
Definition: auxiliary.h:138
int p_LastVblock(poly p, const ring r)
Definition: shiftop.cc:397
#define StringAppend
Definition: emacs.cc:79
int i
Definition: cfEzgcd.cc:125
void PrintS(const char *s)
Definition: reporter.cc:284
BOOLEAN p_LPDivisibleBy(poly a, poly b, const ring r)
Definition: shiftop.cc:688
#define p_LmCheckPolyRing1(p, r)
Definition: monomials.h:177
static unsigned pLength(poly a)
Definition: p_polys.h:191
#define IDELEMS(i)
Definition: simpleideals.h:23
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff &#39;n&#39; represents the zero element.
Definition: coeffs.h:464
#define p_LmTest(p, r)
Definition: p_polys.h:163
#define p_Test(p, r)
Definition: p_polys.h:162
poly p_LPVarAt(poly p, int pos, const ring r)
Definition: shiftop.cc:743
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:856
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent : VarOffset encodes the position in p->exp
Definition: p_polys.h:487
#define pIfThen1(cond, check)
Definition: monomials.h:179
#define NULL
Definition: omList.c:12
b *CanonicalForm B
Definition: facBivar.cc:52
#define R
Definition: sirandom.c:26
const CanonicalForm & w
Definition: facAbsFact.cc:55
int p_FirstVblock(poly p, const ring r)
Definition: shiftop.cc:447
#define pNext(p)
Definition: monomials.h:36
#define p_MemCopy_LengthGeneral(d, s, length)
Definition: p_MemCopy.h:79
#define pSetCoeff0(p, n)
Definition: monomials.h:59
static poly p_GetExp_k_n(poly p, int l, int k, const ring r)
Definition: p_polys.h:1307
void p_mLPunshift(poly m, const ring ri)
Definition: shiftop.cc:322
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1042
poly shift_pp_Mult_mm(poly p, const poly m, const ring ri)
Definition: shiftop.cc:21
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:455
void p_wrp(poly p, ring lmRing, ring tailRing)
Definition: polys0.cc:280
int p
Definition: cfModGcd.cc:4019
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:891
void p_mLPshift(poly m, int sh, const ring ri)
Definition: shiftop.cc:357
ring freeAlgebra(ring r, int d)
create the letterplace ring corresponding to r up to degree d
Definition: shiftop.cc:775
void p_LPunshift(poly p, const ring ri)
Definition: shiftop.cc:348
int BOOLEAN
Definition: auxiliary.h:85
#define pAssume1(cond)
Definition: monomials.h:171
#define omAlloc0(size)
Definition: omAllocDecl.h:211
void k_SplitFrame(poly &m1, poly &m2, int at, const ring r)
Definition: shiftop.cc:585
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
Definition: p_polys.h:291
#define omStrDup(s)
Definition: omAllocDecl.h:263