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 /*
12  * NEEDED BY
13  * - ncHilb.lib
14  */
15 #define SHIFT_MULT_COMPAT_MODE
16 
17 #ifdef SHIFT_MULT_DEBUG
18 #include "../kernel/polys.h"
19 #endif
20 
shift_pp_Mult_mm(poly p,const poly m,const ring ri)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     pNext(q)=NULL;
64     pSetCoeff0(q, n_Mult(mCoeff, pGetCoeff(p), ri->cf));
65 
66     p_GetExpV(p, pExpV, ri);
67     p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
68     p_MemCopy_LengthGeneral(q->exp, p->exp, ri->ExpL_Size); // otherwise q is not initialized correctly
69     p_SetExpV(q, pExpV, ri);
70 
71     pIter(p);
72   }
73   while (p != NULL);
74   omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
75   omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
76   pNext(q) = NULL;
77 #ifdef SHIFT_MULT_COMPAT_MODE
78   p_Delete(&_m, ri); // in this case we copied _m before
79   p_Delete(&pCopyHead, ri); // in this case we copied p before
80 #endif
81 #ifdef SHIFT_MULT_DEBUG
82   PrintLn(); PrintS("shift_pp_Mult_mm result: "); p_wrp(pNext(&rp), ri, ri); PrintLn();
83 #endif
84   p_Test(pNext(&rp), ri);
85   return pNext(&rp);
86 }
87 
88 // destroys p
shift_p_Mult_mm(poly p,const poly m,const ring ri)89 poly shift_p_Mult_mm(poly p, const poly m, const ring ri)
90 {
91 #ifdef SHIFT_MULT_DEBUG
92   PrintLn(); PrintS("shift_p_Mult_mm: ("); p_wrp(p, ri, ri); PrintS(") * "); p_wrp(m, ri, ri);
93 #endif
94 
95   p_Test(p, ri);
96   p_LmTest(m, ri);
97   pAssume(m != NULL);
98   assume(p!=NULL);
99 
100   int lV = ri->isLPring;
101   poly _m = m; // temp hack because m is const
102 #ifdef SHIFT_MULT_COMPAT_MODE
103   _m = p_Copy(_m, ri);
104   p_mLPunshift(_m, ri);
105   p_LPunshift(p, ri);
106 #else
107   assume(p_mFirstVblock(_m, ri) <= 1);
108   assume(p_FirstVblock(p, ri) <= 1); // TODO check that each block is <=1
109 #endif
110   // at this point _m and p are shifted to 1
111 
112   poly q = p; // we use p for iterating and q for the result
113   number mCoeff = pGetCoeff(_m);
114   number pCoeff;
115   pAssume(!n_IsZero(mCoeff, ri->cf));
116 
117   int *mExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
118   p_GetExpV(_m,mExpV,ri);
119   int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
120   int *pExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
121   while (p != NULL)
122   {
123     pCoeff = pGetCoeff(p);
124     pSetCoeff0(p, n_Mult(mCoeff, pCoeff, ri->cf));
125     n_Delete(&pCoeff, ri->cf); // delete the old coeff
126 
127     p_GetExpV(p,pExpV,ri);
128     p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
129     p_SetExpV(p, pExpV, ri);
130 
131     pIter(p);
132   }
133   omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
134   omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
135 #ifdef SHIFT_MULT_COMPAT_MODE
136   p_Delete(&_m, ri); // in this case we copied _m before
137 #endif
138 #ifdef SHIFT_MULT_DEBUG
139   PrintLn(); PrintS("shift_p_Mult_mm result: "); p_wrp(q, ri, ri); PrintLn();
140 #endif
141   p_Test(q, ri);
142   return q;
143 }
144 
shift_pp_mm_Mult(poly p,const poly m,const ring ri)145 poly shift_pp_mm_Mult(poly p, const poly m, const ring ri)
146 {
147 #ifdef SHIFT_MULT_DEBUG
148   PrintLn(); PrintS("shift_pp_mm_Mult: "); p_wrp(m, ri, ri); PrintS(" * ("); p_wrp(p, ri, ri); PrintS(")");
149 #endif
150 
151   p_Test(p, ri);
152   p_LmTest(m, ri);
153   if (p == NULL)
154   {
155     return NULL;
156   }
157 
158   int lV = ri->isLPring;
159   poly _m = m; // temp hack because m is const
160 #ifdef SHIFT_MULT_COMPAT_MODE
161   _m = p_Copy(_m, ri);
162   p_mLPunshift(_m, ri);
163   p = p_Copy(p, ri);
164   poly pCopyHead = p; // used to delete p later
165   p_LPunshift(p, ri);
166 #else
167   assume(p_mFirstVblock(_m, ri) <= 1);
168   assume(p_FirstVblock(p, ri) <= 1); // TODO check that each block is <=1
169 #endif
170   // at this point _m and p are shifted to 1
171 
172   spolyrec rp;
173   poly q = &rp; // we use p for iterating and q for the result
174   number mCoeff = pGetCoeff(_m);
175   omBin bin = ri->PolyBin;
176   pAssume(!n_IsZero(mCoeff, ri->cf));
177   pAssume1(p_GetComp(m, ri) == 0 || p_MaxComp(p, ri) == 0);
178 
179   int *mExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
180   p_GetExpV(_m,mExpV,ri);
181   int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
182   int *pExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
183   do
184   {
185     p_AllocBin(pNext(q), bin, ri);
186     pIter(q);
187     pNext(q)=NULL;
188     pSetCoeff0(q, n_Mult(mCoeff, pGetCoeff(p), ri->cf));
189 
190     p_GetExpV(p, pExpV, ri);
191     p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
192     p_MemCopy_LengthGeneral(q->exp, p->exp, ri->ExpL_Size); // otherwise q is not initialized correctly
193     p_SetExpV(q, pExpV, ri);
194 
195     pIter(p);
196   }
197   while (p != NULL);
198   omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
199   omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
200   pNext(q) = NULL;
201 #ifdef SHIFT_MULT_COMPAT_MODE
202   p_Delete(&_m, ri); // in this case we copied _m before
203   p_Delete(&pCopyHead, ri); // in this case we copied p before
204 #endif
205 #ifdef SHIFT_MULT_DEBUG
206   PrintLn(); PrintS("shift_pp_mm_Mult result: "); p_wrp(pNext(&rp), ri, ri); PrintLn();
207 #endif
208   p_Test(pNext(&rp), ri);
209   return pNext(&rp);
210 }
211 
212 // destroys p
shift_p_mm_Mult(poly p,const poly m,const ring ri)213 poly shift_p_mm_Mult(poly p, const poly m, const ring ri)
214 {
215 #ifdef SHIFT_MULT_DEBUG
216   PrintLn(); PrintS("shift_p_mm_Mult: "); p_wrp(m, ri, ri); PrintS(" * ("); p_wrp(p, ri, ri); PrintS(")");
217 #endif
218 
219   p_Test(p, ri);
220   p_LmTest(m, ri);
221   pAssume(m != NULL);
222   assume(p!=NULL);
223 
224   int lV = ri->isLPring;
225   poly _m = m; // temp hack because m is const
226 #ifdef SHIFT_MULT_COMPAT_MODE
227   _m = p_Copy(_m, ri);
228   p_mLPunshift(_m, ri);
229   p_LPunshift(p, ri);
230 #else
231   assume(p_mFirstVblock(_m, ri) <= 1);
232   assume(p_FirstVblock(p, ri) <= 1); // TODO check that each block is <=1
233 #endif
234   // at this point _m and p are shifted to 1
235 
236   poly q = p; // we use p for iterating and q for the result
237   number mCoeff = pGetCoeff(_m);
238   number pCoeff;
239   pAssume(!n_IsZero(mCoeff, ri->cf));
240 
241   int *mExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
242   p_GetExpV(_m,mExpV,ri);
243   int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
244   int *pExpV = (int *) omAlloc((ri->N+1)*sizeof(int));
245   while (p != NULL)
246   {
247     pCoeff = pGetCoeff(p);
248     pSetCoeff0(p, n_Mult(mCoeff, pCoeff, ri->cf));
249     n_Delete(&pCoeff, ri->cf); // delete the old coeff
250 
251     p_GetExpV(p,pExpV,ri);
252     p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
253     p_SetExpV(p, pExpV, ri);
254 
255     pIter(p);
256   }
257   omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
258   omFreeSize((ADDRESS) mExpV, (ri->N+1)*sizeof(int));
259 #ifdef SHIFT_MULT_COMPAT_MODE
260   p_Delete(&_m, ri); // in this case we copied _m before
261 #endif
262 #ifdef SHIFT_MULT_DEBUG
263   PrintLn(); PrintS("shift_p_mm_Mult result: "); p_wrp(q, ri, ri); PrintLn();
264 #endif
265   p_Test(q, ri);
266   return q;
267 }
268 
269 // p - m*q destroys p
shift_p_Minus_mm_Mult_qq(poly p,poly m,poly q,int & Shorter,const poly spNoether,const ring ri)270 poly shift_p_Minus_mm_Mult_qq(poly p, poly m, poly q, int& Shorter, const poly spNoether, const ring ri) {
271 #ifdef SHIFT_MULT_DEBUG
272   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);
273 #endif
274 
275   Shorter = pLength(p) + pLength(q);
276 
277   poly qq = p_Add_q(p, shift_pp_mm_Mult(q, p_Neg(p_Copy(m, ri), ri), ri), ri);
278 
279 #ifdef SHIFT_MULT_DEBUG
280   PrintLn(); PrintS("shift_p_Minus_mm_Mult_qq result: "); p_wrp(qq, ri, ri); PrintLn();
281 #endif
282   Shorter -= pLength(qq);
283   return qq;
284 }
285 
286 // Unsupported Operation STUBs
shift_pp_Mult_mm_Noether_STUB(poly p,const poly m,const poly spNoether,int & ll,const ring ri)287 poly shift_pp_Mult_mm_Noether_STUB(poly p, const poly m, const poly spNoether, int &ll, const ring ri) {
288   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.");
289 
290   int pLen = 0;
291   if (ll >= 0)
292   {
293     pLen = pLength(p);
294   }
295 
296   p = shift_pp_Mult_mm(p, m, ri);
297 
298   if (ll >= 0)
299   {
300     ll = pLen - pLength(p);
301   }
302   else
303   {
304     ll = pLength(p);
305   }
306 
307   return p;
308 }
309 
310 
shift_pp_Mult_Coeff_mm_DivSelectMult_STUB(poly p,const poly m,const poly a,const poly b,int & shorter,const ring r)311 poly shift_pp_Mult_Coeff_mm_DivSelectMult_STUB(poly p,const poly m, const poly a, const poly b, int &shorter,const ring r) {
312   PrintLn(); WarnS("pp_Mult_Coeff_mm_DivSelectMult is not supported yet by Letterplace. This might lead to unexpected behavior.");
313   return NULL;
314 }
315 
shift_pp_Mult_Coeff_mm_DivSelect_STUB(poly p,const poly m,int & shorter,const ring r)316 poly shift_pp_Mult_Coeff_mm_DivSelect_STUB(poly p, const poly m, int &shorter, const ring r) {
317   PrintLn(); WarnS("pp_Mult_Coeff_mm_DivSelect is not supported yet by Letterplace. This might lead to unexpected behavior.");
318   return NULL;
319 }
320 
321 // auxiliary
322 
323 // unshifts the monomial m
p_mLPunshift(poly m,const ring ri)324 void p_mLPunshift(poly m, const ring ri)
325 {
326   if (m == NULL || p_LmIsConstantComp(m,ri)) return;
327 
328   int lV = ri->isLPring;
329 
330   int shift = p_mFirstVblock(m, ri) - 1;
331 
332   if (shift == 0) return;
333 
334   int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
335   int *s=(int *)omAlloc0((ri->N+1)*sizeof(int));
336   p_GetExpV(m, e, ri);
337 
338   int expVoffset = shift*lV;
339   for (int i = 1 + expVoffset; i <= ri->N; i++)
340   {
341     assume(e[i] <= 1);
342     s[i - expVoffset] = e[i];
343   }
344   p_SetExpV(m,s,ri);
345   omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
346   omFreeSize((ADDRESS) s, (ri->N+1)*sizeof(int));
347 }
348 
349 // unshifts the polynomial p, note: the ordering can be destroyed if the shifts for the monomials are not equal
p_LPunshift(poly p,const ring ri)350 void p_LPunshift(poly p, const ring ri)
351 {
352   while (p!=NULL)
353   {
354     p_mLPunshift(p, ri);
355     pIter(p);
356   }
357 }
358 
p_mLPshift(poly m,int sh,const ring ri)359 void p_mLPshift(poly m, int sh, const ring ri)
360 {
361   if (sh == 0 || m == NULL || p_LmIsConstantComp(m,ri)) return;
362 
363   int lV = ri->isLPring;
364 
365   assume(p_mFirstVblock(m,ri) + sh >= 1);
366   assume(p_mLastVblock(m,ri) + sh <= ri->N/lV);
367 
368   int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
369   int *s=(int *)omAlloc0((ri->N+1)*sizeof(int));
370   p_GetExpV(m,e,ri);
371 
372   if (p_mLastVblock(m, e, ri) + sh > ri->N/lV)
373   {
374     Werror("degree bound of Letterplace ring is %d, but at least %d is needed for this shift", ri->N/lV, p_mLastVblock(m, e, ri) + sh);
375   }
376   for (int i = ri->N - sh*lV; i > 0; i--)
377   {
378     assume(e[i]<=1);
379     if (e[i]==1)
380     {
381       s[i + (sh*lV)] = e[i]; /* actually 1 */
382     }
383   }
384   p_SetExpV(m,s,ri);
385   omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
386   omFreeSize((ADDRESS) s, (ri->N+1)*sizeof(int));
387 }
388 
p_LPshift(poly p,int sh,const ring ri)389 void p_LPshift(poly p, int sh, const ring ri)
390 {
391   if (sh == 0) return;
392 
393   while (p!=NULL)
394   {
395     p_mLPshift(p, sh, ri);
396     pIter(p);
397   }
398 }
399 
400 /* returns the number of maximal block */
401 /* appearing among the monomials of p */
402 /* the 0th block is the 1st one */
p_LastVblock(poly p,const ring r)403 int p_LastVblock(poly p, const ring r)
404 {
405   poly q = p;
406   int ans = 0;
407   while (q!=NULL)
408   {
409     int ansnew = p_mLastVblock(q, r);
410     ans    = si_max(ans,ansnew);
411     pIter(q);
412   }
413   return(ans);
414 }
415 
416 /* for a monomial p, returns the number of the last block */
417 /* where a nonzero exponent is sitting */
p_mLastVblock(poly p,const ring ri)418 int p_mLastVblock(poly p, const ring ri)
419 {
420   if (p == NULL || p_LmIsConstantComp(p,ri))
421   {
422     return(0);
423   }
424 
425   int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
426   p_GetExpV(p,e,ri);
427   int b = p_mLastVblock(p, e, ri);
428   omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
429   return b;
430 }
431 
432 /* for a monomial p with exponent vector expV, returns the number of the last block */
433 /* where a nonzero exponent is sitting */
p_mLastVblock(poly p,int * expV,const ring ri)434 int p_mLastVblock(poly p, int *expV, const ring ri)
435 {
436   if (p == NULL || p_LmIsConstantComp(p,ri))
437   {
438     return(0);
439   }
440 
441   int lV = ri->isLPring;
442   int j,b;
443   j = ri->N;
444   while ( (!expV[j]) && (j>=1) ) j--;
445   assume(j>0);
446   b = (int)((j+lV-1)/lV); /* the number of the block, >=1 */
447   return b;
448 }
449 
450 /* returns the number of maximal block */
451 /* appearing among the monomials of p */
452 /* the 0th block is the 1st one */
p_FirstVblock(poly p,const ring r)453 int p_FirstVblock(poly p, const ring r)
454 {
455   if (p == NULL) {
456     return 0;
457   }
458 
459   poly q = p;
460   int ans = p_mFirstVblock(q, r);
461   while (q!=NULL)
462   {
463     int ansnew = p_mFirstVblock(q, r);
464     if (ansnew > 0) { // don't count constants
465       ans = si_min(ans,ansnew);
466     }
467     pIter(q);
468   }
469   /* do not need to delete q */
470   return(ans);
471 }
472 
473 /* for a monomial p, returns the number of the first block */
474 /* where a nonzero exponent is sitting */
p_mFirstVblock(poly p,const ring ri)475 int p_mFirstVblock(poly p, const ring ri)
476 {
477   if (p == NULL || p_LmIsConstantComp(p,ri))
478   {
479     return(0);
480   }
481 
482   int *e=(int *)omAlloc((ri->N+1)*sizeof(int));
483   p_GetExpV(p,e,ri);
484   int b = p_mFirstVblock(p, e, ri);
485   omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
486   return b;
487 }
488 
489 /* for a monomial p with exponent vector expV, returns the number of the first block */
490 /* where a nonzero exponent is sitting */
p_mFirstVblock(poly p,int * expV,const ring ri)491 int p_mFirstVblock(poly p, int *expV, const ring ri)
492 {
493   if (p == NULL || p_LmIsConstantComp(p,ri))
494   {
495     return(0);
496   }
497 
498   int lV = ri->isLPring;
499   int j,b;
500   j = 1;
501   while ( (!expV[j]) && (j<=ri->N-1) ) j++;
502   assume(j <= ri->N);
503   b = (int)(j+lV-1)/lV; /* the number of the block, 1<= b <= r->N  */
504   return b;
505 }
506 
507 // appends m2ExpV to m1ExpV, also adds their components (one of them is always zero)
p_LPExpVappend(int * m1ExpV,int * m2ExpV,int m1Length,int m2Length,const ring ri)508 void p_LPExpVappend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri) {
509 #ifdef SHIFT_MULT_DEBUG
510   PrintLn(); PrintS("Append");
511   PrintLn(); WriteLPExpV(m1ExpV, ri);
512   PrintLn(); WriteLPExpV(m2ExpV, ri);
513 #endif
514   int last = m1Length + m2Length;
515   if (last > ri->N)
516   {
517     Werror("degree bound of Letterplace ring is %d, but at least %d is needed for this multiplication", ri->N/ri->isLPring, last/ri->isLPring);
518     last = ri->N;
519   }
520   for (int i = 1 + m1Length; i < 1 + last; ++i)
521   {
522     assume(m2ExpV[i - m1Length] <= 1);
523     m1ExpV[i] = m2ExpV[i - m1Length];
524   }
525 
526   assume(m1ExpV[0] == 0 || m2ExpV[0] == 0); // one component should be zero (otherwise this doesn't make any sense)
527   m1ExpV[0] += m2ExpV[0]; // as in the commutative variant (they use MemAdd)
528 #ifdef SHIFT_MULT_DEBUG
529   PrintLn(); WriteLPExpV(m1ExpV, ri);
530 #endif
531   assume(_p_mLPNCGenValid(m1ExpV, ri));
532 }
533 
534 // prepends m2ExpV to m1ExpV, also adds their components (one of them is always zero)
p_LPExpVprepend(int * m1ExpV,int * m2ExpV,int m1Length,int m2Length,const ring ri)535 void p_LPExpVprepend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri)
536 {
537 #ifdef SHIFT_MULT_DEBUG
538   PrintLn(); PrintS("Prepend");
539   PrintLn(); WriteLPExpV(m1ExpV, ri);
540   PrintLn(); WriteLPExpV(m2ExpV, ri);
541 #endif
542   int last = m1Length + m2Length;
543   if (last > ri->N)
544   {
545     Werror("degree bound of Letterplace ring is %d, but at least %d is needed for this multiplication", ri->N/ri->isLPring, last/ri->isLPring);
546     last = ri->N;
547   }
548 
549   // shift m1 by m2Length
550   for (int i = last; i >= 1 + m2Length; --i)
551   {
552     m1ExpV[i] = m1ExpV[i - m2Length];
553   }
554 
555   // write m2 to m1
556   for (int i = 1; i < 1 + m2Length; ++i)
557   {
558     assume(m2ExpV[i] <= 1);
559     m1ExpV[i] = m2ExpV[i];
560   }
561 
562   assume(m1ExpV[0] == 0 || m2ExpV[0] == 0); // one component should be zero (otherwise this doesn't make any sense)
563   m1ExpV[0] += m2ExpV[0]; // as in the commutative variant (they use MemAdd)
564 #ifdef SHIFT_MULT_DEBUG
565   PrintLn(); WriteLPExpV(m1ExpV, ri);
566 #endif
567   assume(_p_mLPNCGenValid(m1ExpV, ri));
568 }
569 
WriteLPExpV(int * expV,ring ri)570 void WriteLPExpV(int *expV, ring ri)
571 {
572   char *s = LPExpVString(expV, ri);
573   PrintS(s);
574   omFree(s);
575 }
576 
LPExpVString(int * expV,ring ri)577 char* LPExpVString(int *expV, ring ri)
578 {
579   StringSetS("");
580   for (int i = 0; i <= ri->N; ++i)
581   {
582     StringAppend("%d", expV[i]);
583     if (i == 0)
584     {
585       StringAppendS("| ");
586     }
587     if (i % ri->isLPring == 0 && i != ri->N)
588     {
589       StringAppendS(" ");
590     }
591   }
592   return StringEndS();
593 }
594 
595 // splits a frame (e.g. x(1)*y(5)) m1 into m1 and m2 (e.g. m1=x(1) and m2=y(1))
596 // at is the number of the block to split at, starting at 1
k_SplitFrame(poly & m1,poly & m2,int at,const ring r)597 void k_SplitFrame(poly &m1, poly &m2, int at, const ring r)
598 {
599   assume(at >= 1);
600   assume(at <= r->N/r->isLPring);
601   int lV = r->isLPring;
602   int split = (lV * (at - 1));
603 
604   m2 = p_GetExp_k_n(m1, 1, split, r);
605   p_SetComp(m2, 0, r); // important, otherwise both m1 and m2 have a component set, this leads to problems later
606   p_Setm(m2, r); // p_mLPunshift also implicitly calls p_Setm(), but just for the case this changes in future.
607   p_mLPunshift(m2, r);
608 
609   m1 = p_Head(m1, r);
610   for(int i = split + 1; i <= r->N; i++)
611   {
612     p_SetExp(m1, i, 0, r);
613   }
614   p_Setm(m1, r);
615 
616   assume(p_FirstVblock(m1,r) <= 1);
617   assume(p_FirstVblock(m2,r) <= 1);
618 }
619 
_p_mLPNCGenValid(poly p,const ring r)620 BOOLEAN _p_mLPNCGenValid(poly p, const ring r)
621 {
622   if (p == NULL) return TRUE;
623   int *e=(int *)omAlloc((r->N+1)*sizeof(int));
624   p_GetExpV(p,e,r);
625   int b = _p_mLPNCGenValid(e, r);
626   omFreeSize((ADDRESS) e, (r->N+1)*sizeof(int));
627   return b;
628 }
629 
_p_mLPNCGenValid(int * mExpV,const ring r)630 BOOLEAN _p_mLPNCGenValid(int *mExpV, const ring r)
631 {
632   BOOLEAN hasNCGen = FALSE;
633   int lV = r->isLPring;
634   int degbound = r->N/lV;
635   int ncGenCount = r->LPncGenCount;
636   for (int i = 1; i <= degbound; i++)
637   {
638     for (int j = i*lV; j > (i*lV - ncGenCount); j--)
639     {
640       if (mExpV[j])
641       {
642         if (hasNCGen)
643         {
644           return FALSE;
645         }
646         hasNCGen = TRUE;
647       }
648     }
649   }
650   return TRUE;
651 }
652 
p_GetNCGen(poly p,const ring r)653 int p_GetNCGen(poly p, const ring r)
654 {
655   if (p == NULL) return 0;
656   assume(_p_mLPNCGenValid(p, r));
657 
658   int lV = r->isLPring;
659   int degbound = r->N/lV;
660   int ncGenCount = r->LPncGenCount;
661   for (int i = 1; i <= degbound; i++)
662   {
663     for (int j = i*lV; j > (i*lV - ncGenCount); j--)
664     {
665       if (p_GetExp(p, j, r))
666       {
667         return j - i*lV + ncGenCount;
668       }
669     }
670   }
671   return 0;
672 }
673 
674 /* tests whether each polynomial of an ideal I lies in in V */
id_IsInV(ideal I,const ring r)675 int id_IsInV(ideal I, const ring r)
676 {
677   int i;
678   int s    = IDELEMS(I)-1;
679   for(i = 0; i <= s; i++)
680   {
681     if ( !p_IsInV(I->m[i], r) )
682     {
683       return(0);
684     }
685   }
686   return(1);
687 }
688 
689 /* tests whether the whole polynomial p in in V */
p_IsInV(poly p,const ring r)690 int p_IsInV(poly p, const ring r)
691 {
692   poly q = p;
693   while (q!=NULL)
694   {
695     if ( !p_mIsInV(q, r) )
696     {
697       return(0);
698     }
699     q = pNext(q);
700   }
701   return(1);
702 }
703 
704 /* there should be two routines: */
705 /* 1. test place-squarefreeness: in homog this suffices: isInV */
706 /* 2. test the presence of a hole -> in the tail??? */
707 
p_mIsInV(poly p,const ring r)708 int p_mIsInV(poly p, const ring r)
709 {
710   int lV = r->isLPring;
711   /* investigate only the leading monomial of p in currRing */
712   if ( p_Totaldegree(p, r)==0 ) return(1);
713   /* returns 1 iff p is in V */
714   /* that is in each block up to a certain one there is only one nonzero exponent */
715   /* lV = the length of V = the number of orig vars */
716   int *e = (int *)omAlloc((r->N+1)*sizeof(int));
717   int  b = (int)((r->N+lV-1)/lV); /* the number of blocks */
718   //int b  = (int)(currRing->N)/lV;
719   int *B = (int *)omAlloc0((b+1)*sizeof(int)); /* the num of elements in a block */
720   p_GetExpV(p,e,r);
721   int i,j;
722   for (j=1; j<=b; j++)
723   {
724     /* we go through all the vars */
725     /* by blocks in lV vars */
726     for (i=(j-1)*lV + 1; i<= j*lV; i++)
727     {
728       if (e[i]) B[j] = B[j]+1;
729     }
730   }
731   //  j = b;
732   //  while ( (!B[j]) && (j>=1)) j--;
733   for (j=b; j>=1; j--)
734   {
735     if (B[j]!=0) break;
736   }
737 
738   if (j==0)
739   {
740     omFreeSize((ADDRESS) e, (r->N+1)*sizeof(int));
741     omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
742     return 1;
743   }
744 
745   if (!_p_mLPNCGenValid(e, r))
746   {
747     omFreeSize((ADDRESS) e, (r->N+1)*sizeof(int));
748     omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
749     return 0;
750   }
751 
752   omFreeSize((ADDRESS) e, (r->N+1)*sizeof(int));
753 
754 //   {
755 //     /* it is a zero exp vector, which is in V */
756 //     freeT(B, b);
757 //     return(1);
758 //   }
759   /* now B[j] != 0 and we test place-squarefreeness */
760   for (; j>=1; j--)
761   {
762     if (B[j]!=1)
763     {
764       omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
765       return 0;
766     }
767   }
768 
769   omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
770   return 1;
771 }
772 
p_LPDivisibleBy(poly a,poly b,const ring r)773 BOOLEAN p_LPDivisibleBy(poly a, poly b, const ring r)
774 {
775   pIfThen1(b!=NULL, p_LmCheckPolyRing1(b, r));
776   pIfThen1(a!=NULL, p_LmCheckPolyRing1(a, r));
777 
778   if (b == NULL) return TRUE;
779   if (a != NULL && (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r)))
780       return _p_LPLmDivisibleByNoComp(a,b,r);
781   return FALSE;
782 }
783 
p_LPLmDivisibleBy(poly a,poly b,const ring r)784 BOOLEAN p_LPLmDivisibleBy(poly a, poly b, const ring r)
785 {
786   p_LmCheckPolyRing1(b, r);
787   pIfThen1(a != NULL, p_LmCheckPolyRing1(b, r));
788   if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
789     return _p_LPLmDivisibleByNoComp(a, b, r);
790   return FALSE;
791 }
792 
_p_LPLmDivisibleByNoComp(poly a,poly b,const ring r)793 BOOLEAN _p_LPLmDivisibleByNoComp(poly a, poly b, const ring r)
794 {
795 #ifdef SHIFT_MULT_COMPAT_MODE
796   a = p_Head(a, r);
797   p_mLPunshift(a, r);
798   b = p_Head(b, r);
799   p_mLPunshift(b, r);
800 #endif
801   int aLastVblock = p_mLastVblock(a, r);
802   int bLastVblock = p_mLastVblock(b, r);
803   for (int i = 0; i <= bLastVblock - aLastVblock; i++)
804   {
805     bool divisible = true;
806     for (int j = 1; j <= aLastVblock * r->isLPring; j++)
807     {
808       if (p_GetExp(a, j, r) > p_GetExp(b, j + (i * r->isLPring), r))
809       {
810         divisible = false;
811         break;
812       }
813     }
814     if (divisible) return TRUE;
815   }
816 #ifdef SHIFT_MULT_COMPAT_MODE
817   p_Delete(&a, r);
818   p_Delete(&b, r);
819 #endif
820   return FALSE;
821 }
822 
p_LPDivisibleBy(ideal I,poly p,ring r)823 BOOLEAN p_LPDivisibleBy(ideal I, poly p, ring r)
824 {
825   for(int i = 0; i < IDELEMS(I); i++)
826   {
827     if (p_LPDivisibleBy(I->m[i], p, r))
828     {
829       return TRUE;
830     }
831   }
832   return FALSE;
833 }
834 
p_LPVarAt(poly p,int pos,const ring r)835 poly p_LPVarAt(poly p, int pos, const ring r)
836 {
837   if (p == NULL || pos < 1 || pos > (r->N / r->isLPring)) return NULL;
838   poly v = p_One(r);
839   for (int i = (pos-1) * r->isLPring + 1; i <= pos * r->isLPring; i++) {
840     if (p_GetExp(p, i, r)) {
841       p_SetExp(v, i - (pos-1) * r->isLPring, 1, r);
842       return v;
843     }
844   }
845   return v;
846 }
847 
848 
849 /*
850 * substitute the n-th variable by e in m
851 * does not destroy m
852 */
p_mLPSubst(poly m,int n,poly e,const ring r)853 poly p_mLPSubst(poly m, int n, poly e, const ring r)
854 {
855   assume(p_GetComp(e, r) == 0);
856   if (m == NULL) return NULL;
857 
858   int lV = r->isLPring;
859   int degbound = r->N/lV;
860 
861   poly result = p_One(r);
862   poly remaining = p_Head(m, r);
863   p_SetComp(result, p_GetComp(remaining, r), r);
864   p_SetComp(remaining, 0, r);
865   for (int i = 0; i < degbound; i++)
866   {
867     int var = n + lV*i;
868     if (p_GetExp(remaining, var, r)) {
869       if (e == NULL) {
870         p_Delete(&result, r);
871         result = NULL;
872         break;
873       }
874       int startOfBlock = 1 + lV*i;
875       int endOfBlock = lV*(i+1);
876 
877       poly left = p_GetExp_k_n(remaining, startOfBlock, r->N, r);
878       p_SetCoeff(left, n_Copy(p_GetCoeff(remaining, r), r->cf), r);
879       p_mLPunshift(left, r);
880 
881       poly right = p_GetExp_k_n(remaining, 1, endOfBlock, r);
882       p_Delete(&remaining, r);
883       remaining = right;
884 
885       left = p_Mult_q(left, p_Copy(e, r), r);
886       result = p_Mult_q(result, left, r);
887     }
888   }
889   if (result == NULL) {
890     return NULL;
891   } else {
892     p_mLPunshift(remaining, r);
893     return p_Mult_q(result, remaining, r);
894   }
895 }
896 
897 /*
898 * also see p_Subst()
899 * substitute the n-th variable by e in p
900 * does not destroy p
901 */
p_LPSubst(poly p,int n,poly e,const ring r)902 poly p_LPSubst(poly p, int n, poly e, const ring r)
903 {
904   poly res = NULL;
905   while (p!=NULL)
906   {
907     res = p_Add_q(res, p_mLPSubst(p, n, e, r), r);
908     pIter(p);
909   }
910   return res;
911 }
912 
913 /// substitute weights from orderings a,wp,Wp
914 /// by d copies of it at position p
freeAlgebra_weights(const ring old_ring,ring new_ring,int p,int d)915 static BOOLEAN freeAlgebra_weights(const ring old_ring, ring new_ring, int p, int d)
916 {
917   omFree(new_ring->wvhdl[p]);
918   int *w=(int*)omAlloc(new_ring->N*sizeof(int));
919   for(int b=0;b<d;b++)
920   {
921     for(int i=old_ring->N-1;i>=0;i--)
922     {
923       if (old_ring->wvhdl[p][i]<-0) return TRUE;
924       w[b*old_ring->N+i]=old_ring->wvhdl[p][i];
925     }
926   }
927   new_ring->wvhdl[p]=w;
928   new_ring->block1[p]=new_ring->N;
929   return FALSE;
930 }
931 
freeAlgebra(ring r,int d,int ncGenCount)932 ring freeAlgebra(ring r, int d, int ncGenCount)
933 {
934   if (ncGenCount) r = rCopy0(r);
935   char *varname=(char *)omAlloc(20);
936   for (int i = 1; i <= ncGenCount; i++)
937   {
938     sprintf(varname, "ncgen(%d)", i);
939     ring save = r;
940     r = rPlusVar(r, varname, 0);
941     if (r==NULL)
942     {
943       omFreeSize(varname, 20);
944       return NULL; /* error in rPlusVar*/
945     }
946     rDelete(save);
947   }
948   omFreeSize(varname, 20);
949   ring R=rCopy0(r);
950   int p;
951   if((r->order[0]==ringorder_C)
952   ||(r->order[0]==ringorder_c))
953     p=1;
954   else
955     p=0;
956   // create R->N
957   R->N=r->N*d;
958   R->wanted_maxExp=7; /* Tst/Manual/letterplace_liftstd.tst*/
959   R->isLPring=r->N;
960   R->LPncGenCount=ncGenCount;
961   // create R->order
962   BOOLEAN has_order_a=FALSE;
963   while (r->order[p]==ringorder_a)
964   {
965     if (freeAlgebra_weights(r,R,p,d))
966     {
967       WerrorS("weights must be positive");
968       return NULL;
969     }
970     has_order_a=TRUE;
971     p++;
972   }
973   R->block1[p]=R->N; /* only dp,Dp,wp,Wp; will be discarded for lp*/
974   switch(r->order[p])
975   {
976     case ringorder_dp:
977     case ringorder_Dp:
978       break;
979     case ringorder_wp:
980     case ringorder_Wp:
981       if (freeAlgebra_weights(r,R,p,d))
982       {
983         WerrorS("weights must be positive");
984         return NULL;
985       }
986       break;
987     case ringorder_lp:
988     case ringorder_rp:
989     {
990       if(has_order_a)
991       {
992         WerrorS("ordering (a(..),lp/rp not implemented for Letterplace rings");
993         return NULL;
994       }
995       int ** wvhdl=(int**)omAlloc0((r->N+3)*sizeof(int*));
996       rRingOrder_t* ord=(rRingOrder_t*)omAlloc0((r->N+3)*sizeof(rRingOrder_t));
997       int* blk0=(int*)omAlloc0((r->N+3)*sizeof(int));
998       int* blk1=(int*)omAlloc0((r->N+3)*sizeof(int));
999       omFree(R->wvhdl);  R->wvhdl=wvhdl;
1000       omFree(R->order);  R->order=ord;
1001       omFree(R->block0); R->block0=blk0;
1002       omFree(R->block1); R->block1=blk1;
1003       for(int i=0;i<r->N;i++)
1004       {
1005         ord[i+p]=ringorder_a;
1006         //Print("entry:%d->a\n",i+p);
1007         blk0[i+p]=1;
1008         blk1[i+p]=R->N;
1009         wvhdl[i+p]=(int*)omAlloc0(R->N*sizeof(int));
1010         for(int j=0;j<d;j++)
1011         {
1012           assume(j*r->N+i<R->N);
1013           if (r->order[p]==ringorder_lp)
1014             wvhdl[i+p][j*r->N+i]=1;
1015          else
1016             wvhdl[i+p][(j+1)*r->N-i-1]=1;
1017         }
1018       }
1019       ord[r->N+p]=r->order[p]; /* lp or rp */
1020       //Print("entry:%d->lp\n",r->N+p);
1021       blk0[r->N+p]=1;
1022       blk1[r->N+p]=R->N;
1023       // copy component order
1024       if (p==1) ord[0]=r->order[0];
1025       else if (p==0) ord[r->N+1]=r->order[1];
1026       else
1027       { // should never happen:
1028         WerrorS("ordering not implemented for Letterplace rings");
1029         return NULL;
1030       }
1031       //if (p==1) PrintS("entry:0 ->c/C\n");
1032       //else if (p==0) Print("entry:%d ->c/C\n",r->N+1);
1033       break;
1034     }
1035     default: WerrorS("ordering not implemented for Letterplace rings");
1036       return NULL;
1037   }
1038   // create R->names
1039   char **names=(char**)omAlloc(R->N*sizeof(char*));
1040   for(int b=0;b<d;b++)
1041   {
1042     for(int i=r->N-1;i>=0;i--)
1043       names[b*r->N+i]=omStrDup(r->names[i]);
1044   }
1045   for(int i=r->N-1;i>=0;i--) omFree(R->names[i]);
1046   omFree(R->names);
1047   R->names=names;
1048 
1049   if (ncGenCount) rDelete(r);
1050   rComplete(R,TRUE);
1051   return R;
1052 }
1053 #endif
1054