source: git/libpolys/polys/monomials/p_polys.cc @ 37c7fc1

fieker-DuValspielwiese
Last change on this file since 37c7fc1 was 37c7fc1, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: fixed a memory double destruction bug in p_Monic
  • Property mode set to 100644
File size: 80.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    p_polys.cc
6 *  Purpose: implementation of ring independent poly procedures?
7 *  Author:  obachman (Olaf Bachmann)
8 *  Created: 8/00
9 *  Version: $Id$
10 *******************************************************************/
11
12#include <ctype.h>
13
14
15#include <omalloc/omalloc.h>
16#include <misc/auxiliary.h>
17#include <misc/options.h>
18#include <misc/intvec.h>
19
20#include <coeffs/longrat.h> // ???
21
22#include "weight.h"
23#include "simpleideals.h"
24
25#include "monomials/ring.h"
26#include "monomials/p_polys.h"
27
28// #include <???/ideals.h>
29// #include <???/int64vec.h>
30
31#ifndef NDEBUG
32// #include <???/febase.h>
33#endif
34
35#ifdef HAVE_PLURAL
36#include "nc/nc.h"
37#include "nc/sca.h"
38#endif
39
40#include "coeffrings.h"
41#ifdef HAVE_FACTORY
42#include "clapsing.h"
43#endif
44
45/***************************************************************
46 *
47 * Completing what needs to be set for the monomial
48 *
49 ***************************************************************/
50// this is special for the syz stuff
51static int* _components = NULL;
52static long* _componentsShifted = NULL;
53static int _componentsExternal = 0;
54
55BOOLEAN pSetm_error=0;
56
57#ifndef NDEBUG
58# define MYTEST 0
59#else /* ifndef NDEBUG */
60# define MYTEST 0
61#endif /* ifndef NDEBUG */
62
63void p_Setm_General(poly p, const ring r)
64{
65  p_LmCheckPolyRing(p, r);
66  int pos=0;
67  if (r->typ!=NULL)
68  {
69    loop
70    {
71      long ord=0;
72      sro_ord* o=&(r->typ[pos]);
73      switch(o->ord_typ)
74      {
75        case ro_dp:
76        {
77          int a,e;
78          a=o->data.dp.start;
79          e=o->data.dp.end;
80          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r);
81          p->exp[o->data.dp.place]=ord;
82          break;
83        }
84        case ro_wp_neg:
85          ord=POLY_NEGWEIGHT_OFFSET;
86          // no break;
87        case ro_wp:
88        {
89          int a,e;
90          a=o->data.wp.start;
91          e=o->data.wp.end;
92          int *w=o->data.wp.weights;
93#if 1
94          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r)*w[i-a];
95#else
96          long ai;
97          int ei,wi;
98          for(int i=a;i<=e;i++)
99          {
100             ei=p_GetExp(p,i,r);
101             wi=w[i-a];
102             ai=ei*wi;
103             if (ai/ei!=wi) pSetm_error=TRUE;
104             ord+=ai;
105             if (ord<ai) pSetm_error=TRUE;
106          }
107#endif
108          p->exp[o->data.wp.place]=ord;
109          break;
110        }
111      case ro_wp64:
112        {
113          int64 ord=0;
114          int a,e;
115          a=o->data.wp64.start;
116          e=o->data.wp64.end;
117          int64 *w=o->data.wp64.weights64;
118          int64 ei,wi,ai;
119          for(int i=a;i<=e;i++)
120          {
121            //Print("exp %d w %d \n",p_GetExp(p,i,r),(int)w[i-a]);
122            //ord+=((int64)p_GetExp(p,i,r))*w[i-a];
123            ei=(int64)p_GetExp(p,i,r);
124            wi=w[i-a];
125            ai=ei*wi;
126            if(ei!=0 && ai/ei!=wi)
127            {
128              pSetm_error=TRUE;
129              #if SIZEOF_LONG == 4
130              Print("ai %lld, wi %lld\n",ai,wi);
131              #else
132              Print("ai %ld, wi %ld\n",ai,wi);
133              #endif
134            }
135            ord+=ai;
136            if (ord<ai)
137            {
138              pSetm_error=TRUE;
139              #if SIZEOF_LONG == 4
140              Print("ai %lld, ord %lld\n",ai,ord);
141              #else
142              Print("ai %ld, ord %ld\n",ai,ord);
143              #endif
144            }
145          }
146          int64 mask=(int64)0x7fffffff;
147          long a_0=(long)(ord&mask); //2^31
148          long a_1=(long)(ord >>31 ); /*(ord/(mask+1));*/
149
150          //Print("mask: %x,  ord: %d,  a_0: %d,  a_1: %d\n"
151          //,(int)mask,(int)ord,(int)a_0,(int)a_1);
152                    //Print("mask: %d",mask);
153
154          p->exp[o->data.wp64.place]=a_1;
155          p->exp[o->data.wp64.place+1]=a_0;
156//            if(p_Setm_error) Print("***************************\n
157//                                    ***************************\n
158//                                    **WARNING: overflow error**\n
159//                                    ***************************\n
160//                                    ***************************\n");
161          break;
162        }
163        case ro_cp:
164        {
165          int a,e;
166          a=o->data.cp.start;
167          e=o->data.cp.end;
168          int pl=o->data.cp.place;
169          for(int i=a;i<=e;i++) { p->exp[pl]=p_GetExp(p,i,r); pl++; }
170          break;
171        }
172        case ro_syzcomp:
173        {
174          int c=p_GetComp(p,r);
175          long sc = c;
176          int* Components = (_componentsExternal ? _components :
177                             o->data.syzcomp.Components);
178          long* ShiftedComponents = (_componentsExternal ? _componentsShifted:
179                                     o->data.syzcomp.ShiftedComponents);
180          if (ShiftedComponents != NULL)
181          {
182            assume(Components != NULL);
183            assume(c == 0 || Components[c] != 0);
184            sc = ShiftedComponents[Components[c]];
185            assume(c == 0 || sc != 0);
186          }
187          p->exp[o->data.syzcomp.place]=sc;
188          break;
189        }
190        case ro_syz:
191        {
192          const unsigned long c = p_GetComp(p, r);
193          const short place = o->data.syz.place;
194          const int limit = o->data.syz.limit;
195         
196          if (c > limit)
197            p->exp[place] = o->data.syz.curr_index;
198          else if (c > 0)
199          {
200            assume( (1 <= c) && (c <= limit) );
201            p->exp[place]= o->data.syz.syz_index[c];
202          }
203          else
204          {
205            assume(c == 0);
206            p->exp[place]= 0;
207          }
208          break;
209        }
210        // Prefix for Induced Schreyer ordering
211        case ro_isTemp: // Do nothing?? (to be removed into suffix later on...?)
212        {
213          assume(p != NULL);
214
215#ifndef NDEBUG
216#if MYTEST
217          Print("p_Setm_General: isTemp ord: pos: %d, p: ", pos);  p_DebugPrint(p, r, r, 1);
218#endif
219#endif
220          int c = p_GetComp(p, r);
221
222          assume( c >= 0 );
223
224          // Let's simulate case ro_syz above....
225          // Should accumulate (by Suffix) and be a level indicator
226          const int* const pVarOffset = o->data.isTemp.pVarOffset;
227
228          assume( pVarOffset != NULL );
229
230          // TODO: Can this be done in the suffix???
231          for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
232          {
233            const int vo = pVarOffset[i];
234            if( vo != -1) // TODO: optimize: can be done once!
235            {
236              // Hans! Please don't break it again! p_SetExp(p, ..., r, vo) is correct:
237              p_SetExp(p, p_GetExp(p, i, r), r, vo); // copy put them verbatim
238              // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct:
239              assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim
240            }
241          }
242#ifndef NDEBUG
243          for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
244          {
245            const int vo = pVarOffset[i];
246            if( vo != -1) // TODO: optimize: can be done once!
247            {
248              // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct:
249              assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim
250            }
251          }
252#if MYTEST
253//          if( p->exp[o->data.isTemp.start] > 0 )
254//          {
255//            PrintS("Initial Value: "); p_DebugPrint(p, r, r, 1);
256//          }
257#endif
258#endif
259          break;
260        }
261
262        // Suffix for Induced Schreyer ordering
263        case ro_is:
264        {
265#ifndef NDEBUG
266#if MYTEST
267          Print("p_Setm_General: ro_is ord: pos: %d, p: ", pos);  p_DebugPrint(p, r, r, 1);
268#endif
269#endif
270
271          assume(p != NULL);
272
273          int c = p_GetComp(p, r);
274
275          assume( c >= 0 );
276          const ideal F = o->data.is.F;
277          const int limit = o->data.is.limit;
278
279          if( F != NULL && c > limit )
280          {
281#ifndef NDEBUG
282#if MYTEST
283            Print("p_Setm_General: ro_is : in rSetm: pos: %d, c: %d >  limit: %d\n", c, pos, limit); // p_DebugPrint(p, r, r, 1);
284#endif
285#endif
286
287            c -= limit;
288            assume( c > 0 );
289            c--;
290
291            assume( c < IDELEMS(F) ); // What about others???
292
293            const poly pp = F->m[c]; // get reference monomial!!!
294
295#ifndef NDEBUG
296#if MYTEST
297            Print("Respective F[c - %d: %d] pp: ", limit, c); 
298            p_DebugPrint(pp, r, r, 1);
299#endif
300#endif
301
302
303            assume(pp != NULL);
304            if(pp == NULL) break;
305
306            const int start = o->data.is.start;
307            const int end = o->data.is.end;
308
309            assume(start <= end);
310             
311//          const int limit = o->data.is.limit;
312          assume( limit >= 0 );
313
314//        const int st = o->data.isTemp.start;       
315
316          if( c > limit )
317            p->exp[start] = 1;
318//          else
319//            p->exp[start] = 0;
320
321             
322#ifndef NDEBUG
323            Print("p_Setm_General: is(-Temp-) :: c: %d, limit: %d, [st:%d] ===>>> %ld\n", c, limit, start, p->exp[start]);
324#endif       
325   
326
327            for( int i = start; i <= end; i++) // v[0] may be here...
328              p->exp[i] += pp->exp[i]; // !!!!!!!! ADD corresponding LT(F)
329
330       
331
332             
333#ifndef NDEBUG
334            const int* const pVarOffset = o->data.is.pVarOffset;
335
336            assume( pVarOffset != NULL );
337
338            for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
339            {
340              const int vo = pVarOffset[i];
341              if( vo != -1) // TODO: optimize: can be done once!
342                // Hans! Please don't break it again! p_GetExp(p/pp, r, vo) is correct:
343                assume( p_GetExp(p, r, vo) == (p_GetExp(p, i, r) + p_GetExp(pp, r, vo)) );
344            }
345            // TODO: how to check this for computed values???
346#endif
347          } else
348          {
349            const int* const pVarOffset = o->data.is.pVarOffset;
350
351            // What about v[0] - component: it will be added later by
352            // suffix!!!
353            // TODO: Test it!
354            const int vo = pVarOffset[0];
355            if( vo != -1 )
356              p->exp[vo] = c; // initial component v[0]!
357
358#ifndef NDEBUG
359#if MYTEST
360            Print("p_Setm_General: ro_is :: c: %d <= limit: %d, vo: %d, exp: %d\n", c, limit, vo, p->exp[vo]);
361            p_DebugPrint(p, r, r, 1);
362#endif       
363#endif       
364          }
365           
366
367          break;
368        }
369        default:
370          dReportError("wrong ord in rSetm:%d\n",o->ord_typ);
371          return;
372      }
373      pos++;
374      if (pos == r->OrdSize) return;
375    }
376  }
377}
378
379void p_Setm_Syz(poly p, ring r, int* Components, long* ShiftedComponents)
380{
381  _components = Components;
382  _componentsShifted = ShiftedComponents;
383  _componentsExternal = 1;
384  p_Setm_General(p, r);
385  _componentsExternal = 0;
386}
387
388// dummy for lp, ls, etc
389void p_Setm_Dummy(poly p, const ring r)
390{
391  p_LmCheckPolyRing(p, r);
392}
393
394// for dp, Dp, ds, etc
395void p_Setm_TotalDegree(poly p, const ring r)
396{
397  p_LmCheckPolyRing(p, r);
398  p->exp[r->pOrdIndex] = p_Totaldegree(p, r);
399}
400
401// for wp, Wp, ws, etc
402void p_Setm_WFirstTotalDegree(poly p, const ring r)
403{
404  p_LmCheckPolyRing(p, r);
405  p->exp[r->pOrdIndex] = p_WFirstTotalDegree(p, r);
406}
407
408p_SetmProc p_GetSetmProc(ring r)
409{
410  // covers lp, rp, ls,
411  if (r->typ == NULL) return p_Setm_Dummy;
412
413  if (r->OrdSize == 1)
414  {
415    if (r->typ[0].ord_typ == ro_dp &&
416        r->typ[0].data.dp.start == 1 &&
417        r->typ[0].data.dp.end == r->N &&
418        r->typ[0].data.dp.place == r->pOrdIndex)
419      return p_Setm_TotalDegree;
420    if (r->typ[0].ord_typ == ro_wp &&
421        r->typ[0].data.wp.start == 1 &&
422        r->typ[0].data.wp.end == r->N &&
423        r->typ[0].data.wp.place == r->pOrdIndex &&
424        r->typ[0].data.wp.weights == r->firstwv)
425      return p_Setm_WFirstTotalDegree;
426  }
427  return p_Setm_General;
428}
429
430
431/* -------------------------------------------------------------------*/
432/* several possibilities for pFDeg: the degree of the head term       */
433
434/* comptible with ordering */
435long p_Deg(poly a, const ring r)
436{
437  p_LmCheckPolyRing(a, r);
438  assume(p_GetOrder(a, r) == p_WTotaldegree(a, r));
439  return p_GetOrder(a, r);
440}
441
442// p_WTotalDegree for weighted orderings
443// whose first block covers all variables
444long p_WFirstTotalDegree(poly p, const ring r)
445{
446  int i;
447  long sum = 0;
448
449  for (i=1; i<= r->firstBlockEnds; i++)
450  {
451    sum += p_GetExp(p, i, r)*r->firstwv[i-1];
452  }
453  return sum;
454}
455
456/*2
457* compute the degree of the leading monomial of p
458* with respect to weigths from the ordering
459* the ordering is not compatible with degree so do not use p->Order
460*/
461long p_WTotaldegree(poly p, const ring r)
462{
463  p_LmCheckPolyRing(p, r);
464  int i, k;
465  long j =0;
466
467  // iterate through each block:
468  for (i=0;r->order[i]!=0;i++)
469  {
470    int b0=r->block0[i];
471    int b1=r->block1[i];
472    switch(r->order[i])
473    {
474      case ringorder_M:
475        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
476        { // in jedem block:
477          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/]*r->OrdSgn;
478        }
479        break;
480      case ringorder_wp:
481      case ringorder_ws:
482      case ringorder_Wp:
483      case ringorder_Ws:
484        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
485        { // in jedem block:
486          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/];
487        }
488        break;
489      case ringorder_lp:
490      case ringorder_ls:
491      case ringorder_rs:
492      case ringorder_dp:
493      case ringorder_ds:
494      case ringorder_Dp:
495      case ringorder_Ds:
496      case ringorder_rp:
497        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
498        {
499          j+= p_GetExp(p,k,r);
500        }
501        break;
502      case ringorder_a64:
503        {
504          int64* w=(int64*)r->wvhdl[i];
505          for (k=0;k<=(b1 /*r->block1[i]*/ - b0 /*r->block0[i]*/);k++)
506          {
507            //there should be added a line which checks if w[k]>2^31
508            j+= p_GetExp(p,k+1, r)*(long)w[k];
509          }
510          //break;
511          return j;
512        }
513      case ringorder_c:
514      case ringorder_C:
515      case ringorder_S:
516      case ringorder_s:
517      case ringorder_IS:
518      case ringorder_aa:
519        break;
520      case ringorder_a:
521        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
522        { // only one line
523          j+= p_GetExp(p,k, r)*r->wvhdl[i][ k- b0 /*r->block0[i]*/];
524        }
525        //break;
526        return j;
527
528#ifndef NDEBUG
529      default:
530        Print("missing order %d in p_WTotaldegree\n",r->order[i]);
531        break;
532#endif
533    }
534  }
535  return  j;
536}
537
538long p_DegW(poly p, const short *w, const ring R)
539{
540  long r=~0L;
541
542  while (p!=NULL)
543  {
544    long t=totaldegreeWecart_IV(p,R,w);
545    if (t>r) r=t;
546    pIter(p);
547  }
548  return r;
549}
550
551int p_Weight(int i, const ring r)
552{
553  if ((r->firstwv==NULL) || (i>r->firstBlockEnds))
554  {
555    return 1;
556  }
557  return r->firstwv[i-1];
558}
559
560long p_WDegree(poly p, const ring r)
561{
562  if (r->firstwv==NULL) return p_Totaldegree(p, r);
563  p_LmCheckPolyRing(p, r);
564  int i;
565  long j =0;
566
567  for(i=1;i<=r->firstBlockEnds;i++)
568    j+=p_GetExp(p, i, r)*r->firstwv[i-1];
569
570  for (;i<=r->N;i++)
571    j+=p_GetExp(p,i, r)*p_Weight(i, r);
572
573  return j;
574}
575
576
577/* ---------------------------------------------------------------------*/
578/* several possibilities for pLDeg: the maximal degree of a monomial in p*/
579/*  compute in l also the pLength of p                                   */
580
581/*2
582* compute the length of a polynomial (in l)
583* and the degree of the monomial with maximal degree: the last one
584*/
585long pLDeg0(poly p,int *l, const ring r)
586{
587  p_CheckPolyRing(p, r);
588  long k= p_GetComp(p, r);
589  int ll=1;
590
591  if (k > 0)
592  {
593    while ((pNext(p)!=NULL) && (p_GetComp(pNext(p), r)==k))
594    {
595      pIter(p);
596      ll++;
597    }
598  }
599  else
600  {
601     while (pNext(p)!=NULL)
602     {
603       pIter(p);
604       ll++;
605     }
606  }
607  *l=ll;
608  return r->pFDeg(p, r);
609}
610
611/*2
612* compute the length of a polynomial (in l)
613* and the degree of the monomial with maximal degree: the last one
614* but search in all components before syzcomp
615*/
616long pLDeg0c(poly p,int *l, const ring r)
617{
618  assume(p!=NULL);
619#ifdef PDEBUG
620  _p_Test(p,r,PDEBUG);
621#endif
622  p_CheckPolyRing(p, r);
623  long o;
624  int ll=1;
625
626  if (! rIsSyzIndexRing(r))
627  {
628    while (pNext(p) != NULL)
629    {
630      pIter(p);
631      ll++;
632    }
633    o = r->pFDeg(p, r);
634  }
635  else
636  {
637    int curr_limit = rGetCurrSyzLimit(r);
638    poly pp = p;
639    while ((p=pNext(p))!=NULL)
640    {
641      if (p_GetComp(p, r)<=curr_limit/*syzComp*/)
642        ll++;
643      else break;
644      pp = p;
645    }
646#ifdef PDEBUG
647    _p_Test(pp,r,PDEBUG);
648#endif
649    o = r->pFDeg(pp, r);
650  }
651  *l=ll;
652  return o;
653}
654
655/*2
656* compute the length of a polynomial (in l)
657* and the degree of the monomial with maximal degree: the first one
658* this works for the polynomial case with degree orderings
659* (both c,dp and dp,c)
660*/
661long pLDegb(poly p,int *l, const ring r)
662{
663  p_CheckPolyRing(p, r);
664  long k= p_GetComp(p, r);
665  long o = r->pFDeg(p, r);
666  int ll=1;
667
668  if (k != 0)
669  {
670    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
671    {
672      ll++;
673    }
674  }
675  else
676  {
677    while ((p=pNext(p)) !=NULL)
678    {
679      ll++;
680    }
681  }
682  *l=ll;
683  return o;
684}
685
686/*2
687* compute the length of a polynomial (in l)
688* and the degree of the monomial with maximal degree:
689* this is NOT the last one, we have to look for it
690*/
691long pLDeg1(poly p,int *l, const ring r)
692{
693  p_CheckPolyRing(p, r);
694  long k= p_GetComp(p, r);
695  int ll=1;
696  long  t,max;
697
698  max=r->pFDeg(p, r);
699  if (k > 0)
700  {
701    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
702    {
703      t=r->pFDeg(p, r);
704      if (t>max) max=t;
705      ll++;
706    }
707  }
708  else
709  {
710    while ((p=pNext(p))!=NULL)
711    {
712      t=r->pFDeg(p, r);
713      if (t>max) max=t;
714      ll++;
715    }
716  }
717  *l=ll;
718  return max;
719}
720
721/*2
722* compute the length of a polynomial (in l)
723* and the degree of the monomial with maximal degree:
724* this is NOT the last one, we have to look for it
725* in all components
726*/
727long pLDeg1c(poly p,int *l, const ring r)
728{
729  p_CheckPolyRing(p, r);
730  int ll=1;
731  long  t,max;
732
733  max=r->pFDeg(p, r);
734  if (rIsSyzIndexRing(r))
735  {
736    long limit = rGetCurrSyzLimit(r);
737    while ((p=pNext(p))!=NULL)
738    {
739      if (p_GetComp(p, r)<=limit)
740      {
741        if ((t=r->pFDeg(p, r))>max) max=t;
742        ll++;
743      }
744      else break;
745    }
746  }
747  else
748  {
749    while ((p=pNext(p))!=NULL)
750    {
751      if ((t=r->pFDeg(p, r))>max) max=t;
752      ll++;
753    }
754  }
755  *l=ll;
756  return max;
757}
758
759// like pLDeg1, only pFDeg == pDeg
760long pLDeg1_Deg(poly p,int *l, const ring r)
761{
762  assume(r->pFDeg == p_Deg);
763  p_CheckPolyRing(p, r);
764  long k= p_GetComp(p, r);
765  int ll=1;
766  long  t,max;
767
768  max=p_GetOrder(p, r);
769  if (k > 0)
770  {
771    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
772    {
773      t=p_GetOrder(p, r);
774      if (t>max) max=t;
775      ll++;
776    }
777  }
778  else
779  {
780    while ((p=pNext(p))!=NULL)
781    {
782      t=p_GetOrder(p, r);
783      if (t>max) max=t;
784      ll++;
785    }
786  }
787  *l=ll;
788  return max;
789}
790
791long pLDeg1c_Deg(poly p,int *l, const ring r)
792{
793  assume(r->pFDeg == p_Deg);
794  p_CheckPolyRing(p, r);
795  int ll=1;
796  long  t,max;
797
798  max=p_GetOrder(p, r);
799  if (rIsSyzIndexRing(r))
800  {
801    long limit = rGetCurrSyzLimit(r);
802    while ((p=pNext(p))!=NULL)
803    {
804      if (p_GetComp(p, r)<=limit)
805      {
806        if ((t=p_GetOrder(p, r))>max) max=t;
807        ll++;
808      }
809      else break;
810    }
811  }
812  else
813  {
814    while ((p=pNext(p))!=NULL)
815    {
816      if ((t=p_GetOrder(p, r))>max) max=t;
817      ll++;
818    }
819  }
820  *l=ll;
821  return max;
822}
823
824// like pLDeg1, only pFDeg == pTotoalDegree
825long pLDeg1_Totaldegree(poly p,int *l, const ring r)
826{
827  p_CheckPolyRing(p, r);
828  long k= p_GetComp(p, r);
829  int ll=1;
830  long  t,max;
831
832  max=p_Totaldegree(p, r);
833  if (k > 0)
834  {
835    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
836    {
837      t=p_Totaldegree(p, r);
838      if (t>max) max=t;
839      ll++;
840    }
841  }
842  else
843  {
844    while ((p=pNext(p))!=NULL)
845    {
846      t=p_Totaldegree(p, r);
847      if (t>max) max=t;
848      ll++;
849    }
850  }
851  *l=ll;
852  return max;
853}
854
855long pLDeg1c_Totaldegree(poly p,int *l, const ring r)
856{
857  p_CheckPolyRing(p, r);
858  int ll=1;
859  long  t,max;
860
861  max=p_Totaldegree(p, r);
862  if (rIsSyzIndexRing(r))
863  {
864    long limit = rGetCurrSyzLimit(r);
865    while ((p=pNext(p))!=NULL)
866    {
867      if (p_GetComp(p, r)<=limit)
868      {
869        if ((t=p_Totaldegree(p, r))>max) max=t;
870        ll++;
871      }
872      else break;
873    }
874  }
875  else
876  {
877    while ((p=pNext(p))!=NULL)
878    {
879      if ((t=p_Totaldegree(p, r))>max) max=t;
880      ll++;
881    }
882  }
883  *l=ll;
884  return max;
885}
886
887// like pLDeg1, only pFDeg == p_WFirstTotalDegree
888long pLDeg1_WFirstTotalDegree(poly p,int *l, const ring r)
889{
890  p_CheckPolyRing(p, r);
891  long k= p_GetComp(p, r);
892  int ll=1;
893  long  t,max;
894
895  max=p_WFirstTotalDegree(p, r);
896  if (k > 0)
897  {
898    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
899    {
900      t=p_WFirstTotalDegree(p, r);
901      if (t>max) max=t;
902      ll++;
903    }
904  }
905  else
906  {
907    while ((p=pNext(p))!=NULL)
908    {
909      t=p_WFirstTotalDegree(p, r);
910      if (t>max) max=t;
911      ll++;
912    }
913  }
914  *l=ll;
915  return max;
916}
917
918long pLDeg1c_WFirstTotalDegree(poly p,int *l, const ring r)
919{
920  p_CheckPolyRing(p, r);
921  int ll=1;
922  long  t,max;
923
924  max=p_WFirstTotalDegree(p, r);
925  if (rIsSyzIndexRing(r))
926  {
927    long limit = rGetCurrSyzLimit(r);
928    while ((p=pNext(p))!=NULL)
929    {
930      if (p_GetComp(p, r)<=limit)
931      {
932        if ((t=p_Totaldegree(p, r))>max) max=t;
933        ll++;
934      }
935      else break;
936    }
937  }
938  else
939  {
940    while ((p=pNext(p))!=NULL)
941    {
942      if ((t=p_Totaldegree(p, r))>max) max=t;
943      ll++;
944    }
945  }
946  *l=ll;
947  return max;
948}
949
950/***************************************************************
951 *
952 * Maximal Exponent business
953 *
954 ***************************************************************/
955
956static inline unsigned long
957p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r,
958              unsigned long number_of_exp)
959{
960  const unsigned long bitmask = r->bitmask;
961  unsigned long ml1 = l1 & bitmask;
962  unsigned long ml2 = l2 & bitmask;
963  unsigned long max = (ml1 > ml2 ? ml1 : ml2);
964  unsigned long j = number_of_exp - 1;
965
966  if (j > 0)
967  {
968    unsigned long mask = bitmask << r->BitsPerExp;
969    while (1)
970    {
971      ml1 = l1 & mask;
972      ml2 = l2 & mask;
973      max |= ((ml1 > ml2 ? ml1 : ml2) & mask);
974      j--;
975      if (j == 0) break;
976      mask = mask << r->BitsPerExp;
977    }
978  }
979  return max;
980}
981
982static inline unsigned long
983p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r)
984{
985  return p_GetMaxExpL2(l1, l2, r, r->ExpPerLong);
986}
987
988poly p_GetMaxExpP(poly p, const ring r)
989{
990  p_CheckPolyRing(p, r);
991  if (p == NULL) return p_Init(r);
992  poly max = p_LmInit(p, r);
993  pIter(p);
994  if (p == NULL) return max;
995  int i, offset;
996  unsigned long l_p, l_max;
997  unsigned long divmask = r->divmask;
998
999  do
1000  {
1001    offset = r->VarL_Offset[0];
1002    l_p = p->exp[offset];
1003    l_max = max->exp[offset];
1004    // do the divisibility trick to find out whether l has an exponent
1005    if (l_p > l_max ||
1006        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1007      max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
1008
1009    for (i=1; i<r->VarL_Size; i++)
1010    {
1011      offset = r->VarL_Offset[i];
1012      l_p = p->exp[offset];
1013      l_max = max->exp[offset];
1014      // do the divisibility trick to find out whether l has an exponent
1015      if (l_p > l_max ||
1016          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1017        max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
1018    }
1019    pIter(p);
1020  }
1021  while (p != NULL);
1022  return max;
1023}
1024
1025unsigned long p_GetMaxExpL(poly p, const ring r, unsigned long l_max)
1026{
1027  unsigned long l_p, divmask = r->divmask;
1028  int i;
1029
1030  while (p != NULL)
1031  {
1032    l_p = p->exp[r->VarL_Offset[0]];
1033    if (l_p > l_max ||
1034        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1035      l_max = p_GetMaxExpL2(l_max, l_p, r);
1036    for (i=1; i<r->VarL_Size; i++)
1037    {
1038      l_p = p->exp[r->VarL_Offset[i]];
1039      // do the divisibility trick to find out whether l has an exponent
1040      if (l_p > l_max ||
1041          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1042        l_max = p_GetMaxExpL2(l_max, l_p, r);
1043    }
1044    pIter(p);
1045  }
1046  return l_max;
1047}
1048
1049
1050
1051
1052/***************************************************************
1053 *
1054 * Misc things
1055 *
1056 ***************************************************************/
1057// returns TRUE, if all monoms have the same component
1058BOOLEAN p_OneComp(poly p, const ring r)
1059{
1060  if(p!=NULL)
1061  {
1062    long i = p_GetComp(p, r);
1063    while (pNext(p)!=NULL)
1064    {
1065      pIter(p);
1066      if(i != p_GetComp(p, r)) return FALSE;
1067    }
1068  }
1069  return TRUE;
1070}
1071
1072/*2
1073*test if a monomial /head term is a pure power
1074*/
1075int p_IsPurePower(const poly p, const ring r)
1076{
1077  int i,k=0;
1078
1079  for (i=r->N;i;i--)
1080  {
1081    if (p_GetExp(p,i, r)!=0)
1082    {
1083      if(k!=0) return 0;
1084      k=i;
1085    }
1086  }
1087  return k;
1088}
1089
1090/*2
1091*test if a polynomial is univariate
1092* return -1 for constant,
1093* 0 for not univariate,s
1094* i if dep. on var(i)
1095*/
1096int p_IsUnivariate(poly p, const ring r)
1097{
1098  int i,k=-1;
1099
1100  while (p!=NULL)
1101  {
1102    for (i=r->N;i;i--)
1103    {
1104      if (p_GetExp(p,i, r)!=0)
1105      {
1106        if((k!=-1)&&(k!=i)) return 0;
1107        k=i;
1108      }
1109    }
1110    pIter(p);
1111  }
1112  return k;
1113}
1114
1115// set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0
1116int  p_GetVariables(poly p, int * e, const ring r)
1117{
1118  int i;
1119  int n=0;
1120  while(p!=NULL)
1121  {
1122    n=0;
1123    for(i=r->N; i>0; i--)
1124    {
1125      if(e[i]==0)
1126      {
1127        if (p_GetExp(p,i,r)>0)
1128        {
1129          e[i]=1;
1130          n++;
1131        }
1132      }
1133      else
1134        n++;
1135    }
1136    if (n==r->N) break;
1137    pIter(p);
1138  }
1139  return n;
1140}
1141
1142
1143/*2
1144* returns a polynomial representing the integer i
1145*/
1146poly p_ISet(int i, const ring r)
1147{
1148  poly rc = NULL;
1149  if (i!=0)
1150  {
1151    rc = p_Init(r);
1152    pSetCoeff0(rc,n_Init(i,r->cf));
1153    if (n_IsZero(pGetCoeff(rc),r->cf))
1154      p_LmDelete(&rc,r);
1155  }
1156  return rc;
1157}
1158
1159/*2
1160* an optimized version of p_ISet for the special case 1
1161*/
1162poly p_One(const ring r)
1163{
1164  poly rc = p_Init(r);
1165  pSetCoeff0(rc,n_Init(1,r->cf));
1166  return rc;
1167}
1168
1169void p_Split(poly p, poly *h)
1170{
1171  *h=pNext(p);
1172  pNext(p)=NULL;
1173}
1174
1175/*2
1176* pair has no common factor ? or is no polynomial
1177*/
1178BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r)
1179{
1180
1181  if (p_GetComp(p1,r) > 0 || p_GetComp(p2,r) > 0)
1182    return FALSE;
1183  int i = rVar(r);
1184  loop
1185  {
1186    if ((p_GetExp(p1, i, r) > 0) && (p_GetExp(p2, i, r) > 0))
1187      return FALSE;
1188    i--;
1189    if (i == 0)
1190      return TRUE;
1191  }
1192}
1193
1194/*2
1195* convert monomial given as string to poly, e.g. 1x3y5z
1196*/
1197const char * p_Read(const char *st, poly &rc, const ring r)
1198{
1199  if (r==NULL) { rc=NULL;return st;}
1200  int i,j;
1201  rc = p_Init(r);
1202  const char *s = r->cf->cfRead(st,&(rc->coef),r->cf);
1203  if (s==st)
1204  /* i.e. it does not start with a coeff: test if it is a ringvar*/
1205  {
1206    j = r_IsRingVar(s,r);
1207    if (j >= 0)
1208    {
1209      p_IncrExp(rc,1+j,r);
1210      while (*s!='\0') s++;
1211      goto done;
1212    }
1213  }
1214  while (*s!='\0')
1215  {
1216    char ss[2];
1217    ss[0] = *s++;
1218    ss[1] = '\0';
1219    j = r_IsRingVar(ss,r);
1220    if (j >= 0)
1221    {
1222      const char *s_save=s;
1223      s = eati(s,&i);
1224      if (((unsigned long)i) >  r->bitmask)
1225      {
1226        // exponent to large: it is not a monomial
1227        p_LmDelete(&rc,r);
1228        return s_save;
1229      }
1230      p_AddExp(rc,1+j, (long)i, r);
1231    }
1232    else
1233    {
1234      // 1st char of is not a varname
1235      p_LmDelete(&rc,r);
1236      s--;
1237      return s;
1238    }
1239  }
1240done:
1241  if (n_IsZero(pGetCoeff(rc),r->cf)) p_LmDelete(&rc,r);
1242  else
1243  {
1244#ifdef HAVE_PLURAL
1245    // in super-commutative ring
1246    // squares of anti-commutative variables are zeroes!
1247    if(rIsSCA(r))
1248    {
1249      const unsigned int iFirstAltVar = scaFirstAltVar(r);
1250      const unsigned int iLastAltVar  = scaLastAltVar(r);
1251
1252      assume(rc != NULL);
1253
1254      for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
1255        if( p_GetExp(rc, k, r) > 1 )
1256        {
1257          p_LmDelete(&rc, r);
1258          goto finish;
1259        }
1260    }
1261#endif
1262
1263    p_Setm(rc,r);
1264  }
1265finish:
1266  return s;
1267}
1268poly p_mInit(const char *st, BOOLEAN &ok, const ring r)
1269{
1270  poly p;
1271  const char *s=p_Read(st,p,r);
1272  if (*s!='\0')
1273  {
1274    if ((s!=st)&&isdigit(st[0]))
1275    {
1276      errorreported=TRUE;
1277    }
1278    ok=FALSE;
1279    p_Delete(&p,r);
1280    return NULL;
1281  }
1282  #ifdef PDEBUG
1283  _p_Test(p,r,PDEBUG);
1284  #endif
1285  ok=!errorreported;
1286  return p;
1287}
1288
1289/*2
1290* returns a polynomial representing the number n
1291* destroys n
1292*/
1293poly p_NSet(number n, const ring r)
1294{
1295  if (n_IsZero(n,r->cf))
1296  {
1297    n_Delete(&n, r->cf);
1298    return NULL;
1299  }
1300  else
1301  {
1302    poly rc = p_Init(r);
1303    pSetCoeff0(rc,n);
1304    return rc;
1305  }
1306}
1307/*2
1308* assumes that the head term of b is a multiple of the head term of a
1309* and return the multiplicant *m
1310* Frank's observation: If LM(b) = LM(a)*m, then we may actually set
1311* negative(!) exponents in the below loop. I suspect that the correct
1312* comment should be "assumes that LM(a) = LM(b)*m, for some monomial m..."
1313*/
1314poly p_Divide(poly a, poly b, const ring r)
1315{
1316  assume((p_GetComp(a,r)==p_GetComp(b,r)) || (p_GetComp(b,r)==0));
1317  int i;
1318  poly result = p_Init(r);
1319
1320  for(i=(int)r->N; i; i--)
1321    p_SetExp(result,i, p_GetExp(a,i,r)- p_GetExp(b,i,r),r);
1322  p_SetComp(result, p_GetComp(a,r) - p_GetComp(b,r),r);
1323  p_Setm(result,r);
1324  return result;
1325}
1326
1327poly p_Div_nn(poly p, const number n, const ring r)
1328{
1329  pAssume(!n_IsZero(n,r->cf));
1330  p_Test(p, r);
1331
1332  poly q = p;
1333  while (p != NULL)
1334  {
1335    number nc = pGetCoeff(p);
1336    pSetCoeff0(p, n_Div(nc, n, r->cf));
1337    n_Delete(&nc, r->cf);
1338    pIter(p);
1339  }
1340  p_Test(q, r);
1341  return q;
1342}
1343
1344/*2
1345* divides a by the monomial b, ignores monomials which are not divisible
1346* assumes that b is not NULL, destroyes b
1347*/
1348poly p_DivideM(poly a, poly b, const ring r)
1349{
1350  if (a==NULL) { p_Delete(&b,r); return NULL; }
1351  poly result=a;
1352  poly prev=NULL;
1353  int i;
1354#ifdef HAVE_RINGS
1355  number inv=pGetCoeff(b);
1356#else
1357  number inv=n_Invers(pGetCoeff(b),r->cf);
1358#endif
1359
1360  while (a!=NULL)
1361  {
1362    if (p_DivisibleBy(b,a,r))
1363    {
1364      for(i=(int)r->N; i; i--)
1365         p_SubExp(a,i, p_GetExp(b,i,r),r);
1366      p_SubComp(a, p_GetComp(b,r),r);
1367      p_Setm(a,r);
1368      prev=a;
1369      pIter(a);
1370    }
1371    else
1372    {
1373      if (prev==NULL)
1374      {
1375        p_LmDelete(&result,r);
1376        a=result;
1377      }
1378      else
1379      {
1380        p_LmDelete(&pNext(prev),r);
1381        a=pNext(prev);
1382      }
1383    }
1384  }
1385#ifdef HAVE_RINGS
1386  if (n_IsUnit(inv,r->cf))
1387  {
1388    inv = n_Invers(inv,r->cf);
1389    p_Mult_nn(result,inv,r);
1390    n_Delete(&inv, r->cf);
1391  }
1392  else
1393  {
1394    p_Div_nn(result,inv,r);
1395  }
1396#else
1397  p_Mult_nn(result,inv,r);
1398  n_Delete(&inv, r->cf);
1399#endif
1400  p_Delete(&b, r);
1401  return result;
1402}
1403
1404#ifdef HAVE_RINGS
1405/* TRUE iff LT(f) | LT(g) */
1406BOOLEAN p_DivisibleByRingCase(poly f, poly g, const ring r)
1407{
1408  int exponent;
1409  for(int i = (int)rVar(r); i>0; i--)
1410  {
1411    exponent = p_GetExp(g, i, r) - p_GetExp(f, i, r);
1412    if (exponent < 0) return FALSE;
1413  }
1414  return n_DivBy(pGetCoeff(g), pGetCoeff(f), r->cf);
1415}
1416#endif
1417
1418/*2
1419* returns the LCM of the head terms of a and b in *m
1420*/
1421void p_Lcm(poly a, poly b, poly m, const ring r)
1422{
1423  int i;
1424  for (i=rVar(r); i; i--)
1425  {
1426    p_SetExp(m,i, si_max( p_GetExp(a,i,r), p_GetExp(b,i,r)),r);
1427  }
1428  p_SetComp(m, si_max(p_GetComp(a,r), p_GetComp(b,r)),r);
1429  /* Don't do a pSetm here, otherwise hres/lres chockes */
1430}
1431
1432/* assumes that p and divisor are univariate polynomials in r,
1433   mentioning the same variable;
1434   assumes divisor != NULL;
1435   p may be NULL;
1436   assumes a global monomial ordering in r;
1437   performs polynomial division of p by divisor:
1438     - afterwards p contains the remainder of the division, i.e.,
1439       p_before = result * divisor + p_afterwards;
1440     - if needResult == TRUE, then the method computes and returns 'result',
1441       otherwise NULL is returned (This parametrization can be used when
1442       one is only interested in the remainder of the division. In this
1443       case, the method will be slightly faster.)
1444   leaves divisor unmodified */
1445poly p_PolyDiv(poly &p, poly divisor, BOOLEAN needResult, ring r)
1446{
1447//printf("p_PolyDiv:\n");
1448  assume(divisor != NULL);
1449  if (p == NULL) return NULL;
1450 
1451  poly result = NULL;
1452  number divisorLC = p_GetCoeff(divisor, r);
1453  int divisorLE = p_GetExp(divisor, 1, r);
1454//p_Write(p, r); p_Write(divisor, r);
1455  while ((p != NULL) && (p_Deg(p, r) >= p_Deg(divisor, r)))
1456  {
1457    /* determine t = LT(p) / LT(divisor) */
1458    poly t = p_ISet(1, r);
1459    number c = n_Div(p_GetCoeff(p, r), divisorLC, r->cf);
1460    p_SetCoeff(t, c, r);
1461    int e = p_GetExp(p, 1, r) - divisorLE;
1462    p_SetExp(t, 1, e, r);
1463    p_Setm(t, r);
1464//printf("t\n");
1465//p_Write(t, r);
1466    if (needResult) result = p_Add_q(result, p_Copy(t, r), r);
1467    p = p_Add_q(p, p_Neg(p_Mult_q(t, p_Copy(divisor, r), r), r), r);
1468//printf("p\n");
1469//p_Write(p, r);
1470  }
1471  return result;
1472}
1473
1474/* returns NULL if p == NULL, otherwise makes p monic by dividing
1475   by its leading coefficient (only done if this is not already 1);
1476   this assumes that we are over a ground field so that division
1477   is well-defined;
1478   modifies p */
1479void p_Monic(poly &p, ring r)
1480{
1481  if (p == NULL) return;
1482  poly pp = p;
1483  number lc = p_GetCoeff(p, r);
1484  if (n_IsOne(lc, r->cf)) return;
1485
1486  if( pNext(p) == NULL )
1487  { 
1488     p_SetCoeff(p, n_Init(1, r->cf), r); // will delete the leading coeff 1st!
1489     return;
1490  }
1491   
1492  // TODO: consider unsing 1xn_Inverse + n_Mult instead of multiple n_Div
1493  p_SetCoeff0(p, n_Init(1, r->cf), r); // no coeff destruction!
1494   
1495  p = pIter(p);
1496   
1497  while (p != NULL)
1498  {
1499    number c = p_GetCoeff(p, r);
1500    number n = n_Div(c, lc, r->cf);
1501    n_Delete(&c, r->cf);
1502    p_SetCoeff(p, n, r);
1503    p = pIter(p);
1504  }
1505  n_Delete(&lc, r->cf);
1506  p = pp;
1507}
1508
1509/* see p_Gcd;
1510   additional assumption: deg(p) >= deg(q);
1511   must destroy p and q (unless one of them is returned) */
1512poly p_GcdHelper(poly &p, poly &q, ring r)
1513{
1514  if (q == NULL)
1515  {
1516    /* We have to make p monic before we return it, so that if the
1517       gcd is a unit in the ground field, we will actually return 1. */
1518    p_Monic(p, r);
1519    return p;
1520  }
1521  else
1522  {
1523    p_PolyDiv(p, q, FALSE, r);
1524    return p_GcdHelper(q, p, r);
1525  }
1526}
1527
1528/* assumes that p and q are univariate polynomials in r,
1529   mentioning the same variable;
1530   assumes a global monomial ordering in r;
1531   assumes that not both p and q are NULL;
1532   returns the gcd of p and q;
1533   leaves p and q unmodified */
1534poly p_Gcd(poly p, poly q, ring r)
1535{
1536  assume((p != NULL) || (q != NULL));
1537 
1538  poly a = p; poly b = q;
1539  if (p_Deg(a, r) < p_Deg(b, r)) { a = q; b = p; }
1540  a = p_Copy(a, r); b = p_Copy(b, r);
1541  return p_GcdHelper(a, b, r);
1542}
1543
1544/* see p_ExtGcd;
1545   additional assumption: deg(p) >= deg(q);
1546   must destroy p and q (unless one of them is returned) */
1547poly p_ExtGcdHelper(poly &p, poly &pFactor, poly &q, poly &qFactor,
1548                    ring r)
1549{
1550  if (q == NULL)
1551  {
1552    qFactor = NULL;
1553    pFactor = p_ISet(1, r);
1554    number n = p_GetCoeff(pFactor, r);
1555//printf("p_ExtGcdHelper0:\n");
1556//p_Write(p, r);
1557    p_SetCoeff(pFactor, n_Invers(p_GetCoeff(p, r), r->cf), r);
1558    n_Delete(&n, r->cf);
1559    p_Monic(p, r);
1560    return p;
1561  }
1562  else
1563  {
1564//printf("p_ExtGcdHelper1:\n");
1565//p_Write(p, r); p_Write(q, r);
1566    poly pDivQ = p_PolyDiv(p, q, TRUE, r);
1567    poly ppFactor = NULL; poly qqFactor = NULL;
1568//printf("p_ExtGcdHelper2:\n");
1569//p_Write(p, r); p_Write(ppFactor, r);
1570//p_Write(q, r); p_Write(qqFactor, r);
1571    poly theGcd = p_ExtGcdHelper(q, qqFactor, p, ppFactor, r);
1572//printf("p_ExtGcdHelper3:\n");
1573//p_Write(q, r); p_Write(qqFactor, r);
1574//p_Write(p, r); p_Write(ppFactor, r);
1575//p_Write(theGcd, r);
1576    pFactor = ppFactor;
1577    qFactor = p_Add_q(qqFactor,
1578                      p_Neg(p_Mult_q(pDivQ, p_Copy(ppFactor, r), r), r),
1579                      r);
1580//printf("p_ExtGcdHelper4:\n");
1581//p_Write(pFactor, r); p_Write(qFactor, r);
1582    return theGcd;
1583  }
1584}
1585
1586/* assumes that p and q are univariate polynomials in r,
1587   mentioning the same variable;
1588   assumes a global monomial ordering in r;
1589   assumes that not both p and q are NULL;
1590   returns the gcd of p and q;
1591   moreover, afterwards pFactor and qFactor contain appropriate
1592   factors such that gcd(p, q) = p * pFactor + q * qFactor;
1593   leaves p and q unmodified */
1594poly p_ExtGcd(poly p, poly &pFactor, poly q, poly &qFactor, ring r)
1595{
1596//printf("p_ExtGcd1:\n");
1597//p_Write(p, r); p_Write(q, r);
1598  assume((p != NULL) || (q != NULL)); 
1599  poly a = p; poly b = q; BOOLEAN aCorrespondsToP = TRUE;
1600  if (p_Deg(a, r) < p_Deg(b, r))
1601  { a = q; b = p; aCorrespondsToP = FALSE; }
1602  a = p_Copy(a, r); b = p_Copy(b, r);
1603  poly aFactor = NULL; poly bFactor = NULL;
1604  poly theGcd = p_ExtGcdHelper(a, aFactor, b, bFactor, r);
1605  if (aCorrespondsToP) { pFactor = aFactor; qFactor = bFactor; }
1606  else                 { pFactor = bFactor; qFactor = aFactor; }
1607  return theGcd;
1608}
1609
1610/*2
1611* returns the partial differentiate of a by the k-th variable
1612* does not destroy the input
1613*/
1614poly p_Diff(poly a, int k, const ring r)
1615{
1616  poly res, f, last;
1617  number t;
1618
1619  last = res = NULL;
1620  while (a!=NULL)
1621  {
1622    if (p_GetExp(a,k,r)!=0)
1623    {
1624      f = p_LmInit(a,r);
1625      t = n_Init(p_GetExp(a,k,r),r->cf);
1626      pSetCoeff0(f,n_Mult(t,pGetCoeff(a),r->cf));
1627      n_Delete(&t,r->cf);
1628      if (n_IsZero(pGetCoeff(f),r->cf))
1629        p_LmDelete(&f,r);
1630      else
1631      {
1632        p_DecrExp(f,k,r);
1633        p_Setm(f,r);
1634        if (res==NULL)
1635        {
1636          res=last=f;
1637        }
1638        else
1639        {
1640          pNext(last)=f;
1641          last=f;
1642        }
1643      }
1644    }
1645    pIter(a);
1646  }
1647  return res;
1648}
1649
1650static poly p_DiffOpM(poly a, poly b,BOOLEAN multiply, const ring r)
1651{
1652  int i,j,s;
1653  number n,h,hh;
1654  poly p=p_One(r);
1655  n=n_Mult(pGetCoeff(a),pGetCoeff(b),r->cf);
1656  for(i=rVar(r);i>0;i--)
1657  {
1658    s=p_GetExp(b,i,r);
1659    if (s<p_GetExp(a,i,r))
1660    {
1661      n_Delete(&n,r->cf);
1662      p_LmDelete(&p,r);
1663      return NULL;
1664    }
1665    if (multiply)
1666    {
1667      for(j=p_GetExp(a,i,r); j>0;j--)
1668      {
1669        h = n_Init(s,r->cf);
1670        hh=n_Mult(n,h,r->cf);
1671        n_Delete(&h,r->cf);
1672        n_Delete(&n,r->cf);
1673        n=hh;
1674        s--;
1675      }
1676      p_SetExp(p,i,s,r);
1677    }
1678    else
1679    {
1680      p_SetExp(p,i,s-p_GetExp(a,i,r),r);
1681    }
1682  }
1683  p_Setm(p,r);
1684  /*if (multiply)*/ p_SetCoeff(p,n,r);
1685  if (n_IsZero(n,r->cf))  p=p_LmDeleteAndNext(p,r); // return NULL as p is a monomial
1686  return p;
1687}
1688
1689poly p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r)
1690{
1691  poly result=NULL;
1692  poly h;
1693  for(;a!=NULL;pIter(a))
1694  {
1695    for(h=b;h!=NULL;pIter(h))
1696    {
1697      result=p_Add_q(result,p_DiffOpM(a,h,multiply,r),r);
1698    }
1699  }
1700  return result;
1701}
1702/*2
1703* subtract p2 from p1, p1 and p2 are destroyed
1704* do not put attention on speed: the procedure is only used in the interpreter
1705*/
1706poly p_Sub(poly p1, poly p2, const ring r)
1707{
1708  return p_Add_q(p1, p_Neg(p2,r),r);
1709}
1710
1711/*3
1712* compute for a monomial m
1713* the power m^exp, exp > 1
1714* destroys p
1715*/
1716static poly p_MonPower(poly p, int exp, const ring r)
1717{
1718  int i;
1719
1720  if(!n_IsOne(pGetCoeff(p),r->cf))
1721  {
1722    number x, y;
1723    y = pGetCoeff(p);
1724    n_Power(y,exp,&x,r->cf);
1725    n_Delete(&y,r->cf);
1726    pSetCoeff0(p,x);
1727  }
1728  for (i=rVar(r); i!=0; i--)
1729  {
1730    p_MultExp(p,i, exp,r);
1731  }
1732  p_Setm(p,r);
1733  return p;
1734}
1735
1736/*3
1737* compute for monomials p*q
1738* destroys p, keeps q
1739*/
1740static void p_MonMult(poly p, poly q, const ring r)
1741{
1742  number x, y;
1743
1744  y = pGetCoeff(p);
1745  x = n_Mult(y,pGetCoeff(q),r->cf);
1746  n_Delete(&y,r->cf);
1747  pSetCoeff0(p,x);
1748  //for (int i=pVariables; i!=0; i--)
1749  //{
1750  //  pAddExp(p,i, pGetExp(q,i));
1751  //}
1752  //p->Order += q->Order;
1753  p_ExpVectorAdd(p,q,r);
1754}
1755
1756/*3
1757* compute for monomials p*q
1758* keeps p, q
1759*/
1760static poly p_MonMultC(poly p, poly q, const ring rr)
1761{
1762  number x;
1763  poly r = p_Init(rr);
1764
1765  x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr->cf);
1766  pSetCoeff0(r,x);
1767  p_ExpVectorSum(r,p, q, rr);
1768  return r;
1769}
1770
1771/*3
1772*  create binomial coef.
1773*/
1774static number* pnBin(int exp, const ring r)
1775{
1776  int e, i, h;
1777  number x, y, *bin=NULL;
1778
1779  x = n_Init(exp,r->cf);
1780  if (n_IsZero(x,r->cf))
1781  {
1782    n_Delete(&x,r->cf);
1783    return bin;
1784  }
1785  h = (exp >> 1) + 1;
1786  bin = (number *)omAlloc0(h*sizeof(number));
1787  bin[1] = x;
1788  if (exp < 4)
1789    return bin;
1790  i = exp - 1;
1791  for (e=2; e<h; e++)
1792  {
1793      x = n_Init(i,r->cf);
1794      i--;
1795      y = n_Mult(x,bin[e-1],r->cf);
1796      n_Delete(&x,r->cf);
1797      x = n_Init(e,r->cf);
1798      bin[e] = n_IntDiv(y,x,r->cf);
1799      n_Delete(&x,r->cf);
1800      n_Delete(&y,r->cf);
1801  }
1802  return bin;
1803}
1804
1805static void pnFreeBin(number *bin, int exp,const coeffs r)
1806{
1807  int e, h = (exp >> 1) + 1;
1808
1809  if (bin[1] != NULL)
1810  {
1811    for (e=1; e<h; e++)
1812      n_Delete(&(bin[e]),r);
1813  }
1814  omFreeSize((ADDRESS)bin, h*sizeof(number));
1815}
1816
1817/*
1818*  compute for a poly p = head+tail, tail is monomial
1819*          (head + tail)^exp, exp > 1
1820*          with binomial coef.
1821*/
1822static poly p_TwoMonPower(poly p, int exp, const ring r)
1823{
1824  int eh, e;
1825  long al;
1826  poly *a;
1827  poly tail, b, res, h;
1828  number x;
1829  number *bin = pnBin(exp,r);
1830
1831  tail = pNext(p);
1832  if (bin == NULL)
1833  {
1834    p_MonPower(p,exp,r);
1835    p_MonPower(tail,exp,r);
1836#ifdef PDEBUG
1837    p_Test(p,r);
1838#endif
1839    return p;
1840  }
1841  eh = exp >> 1;
1842  al = (exp + 1) * sizeof(poly);
1843  a = (poly *)omAlloc(al);
1844  a[1] = p;
1845  for (e=1; e<exp; e++)
1846  {
1847    a[e+1] = p_MonMultC(a[e],p,r);
1848  }
1849  res = a[exp];
1850  b = p_Head(tail,r);
1851  for (e=exp-1; e>eh; e--)
1852  {
1853    h = a[e];
1854    x = n_Mult(bin[exp-e],pGetCoeff(h),r->cf);
1855    p_SetCoeff(h,x,r);
1856    p_MonMult(h,b,r);
1857    res = pNext(res) = h;
1858    p_MonMult(b,tail,r);
1859  }
1860  for (e=eh; e!=0; e--)
1861  {
1862    h = a[e];
1863    x = n_Mult(bin[e],pGetCoeff(h),r->cf);
1864    p_SetCoeff(h,x,r);
1865    p_MonMult(h,b,r);
1866    res = pNext(res) = h;
1867    p_MonMult(b,tail,r);
1868  }
1869  p_LmDelete(&tail,r);
1870  pNext(res) = b;
1871  pNext(b) = NULL;
1872  res = a[exp];
1873  omFreeSize((ADDRESS)a, al);
1874  pnFreeBin(bin, exp, r->cf);
1875//  tail=res;
1876// while((tail!=NULL)&&(pNext(tail)!=NULL))
1877// {
1878//   if(nIsZero(pGetCoeff(pNext(tail))))
1879//   {
1880//     pLmDelete(&pNext(tail));
1881//   }
1882//   else
1883//     pIter(tail);
1884// }
1885#ifdef PDEBUG
1886  p_Test(res,r);
1887#endif
1888  return res;
1889}
1890
1891static poly p_Pow(poly p, int i, const ring r)
1892{
1893  poly rc = p_Copy(p,r);
1894  i -= 2;
1895  do
1896  {
1897    rc = p_Mult_q(rc,p_Copy(p,r),r);
1898    p_Normalize(rc,r);
1899    i--;
1900  }
1901  while (i != 0);
1902  return p_Mult_q(rc,p,r);
1903}
1904
1905/*2
1906* returns the i-th power of p
1907* p will be destroyed
1908*/
1909poly p_Power(poly p, int i, const ring r)
1910{
1911  poly rc=NULL;
1912
1913  if (i==0)
1914  {
1915    p_Delete(&p,r);
1916    return p_One(r);
1917  }
1918
1919  if(p!=NULL)
1920  {
1921    if ( (i > 0) && ((unsigned long ) i > (r->bitmask)))
1922    {
1923      Werror("exponent %d is too large, max. is %ld",i,r->bitmask);
1924      return NULL;
1925    }
1926    switch (i)
1927    {
1928// cannot happen, see above
1929//      case 0:
1930//      {
1931//        rc=pOne();
1932//        pDelete(&p);
1933//        break;
1934//      }
1935      case 1:
1936        rc=p;
1937        break;
1938      case 2:
1939        rc=p_Mult_q(p_Copy(p,r),p,r);
1940        break;
1941      default:
1942        if (i < 0)
1943        {
1944          p_Delete(&p,r);
1945          return NULL;
1946        }
1947        else
1948        {
1949#ifdef HAVE_PLURAL
1950          if (rIsPluralRing(r)) /* in the NC case nothing helps :-( */
1951          {
1952            int j=i;
1953            rc = p_Copy(p,r);
1954            while (j>1)
1955            {
1956              rc = p_Mult_q(p_Copy(p,r),rc,r);
1957              j--;
1958            }
1959            p_Delete(&p,r);
1960            return rc;
1961          }
1962#endif
1963          rc = pNext(p);
1964          if (rc == NULL)
1965            return p_MonPower(p,i,r);
1966          /* else: binom ?*/
1967          int char_p=rChar(r);
1968          if ((pNext(rc) != NULL)
1969#ifdef HAVE_RINGS
1970             || rField_is_Ring(r)
1971#endif
1972             )
1973            return p_Pow(p,i,r);
1974          if ((char_p==0) || (i<=char_p))
1975            return p_TwoMonPower(p,i,r);
1976          poly p_p=p_TwoMonPower(p_Copy(p,r),char_p,r);
1977          return p_Mult_q(p_Power(p,i-char_p,r),p_p,r);
1978        }
1979      /*end default:*/
1980    }
1981  }
1982  return rc;
1983}
1984
1985/* --------------------------------------------------------------------------------*/
1986/* content suff                                                                   */
1987
1988static number p_InitContent(poly ph, const ring r);
1989static number p_InitContent_a(poly ph, const ring r);
1990
1991void p_Content(poly ph, const ring r)
1992{
1993#ifdef HAVE_RINGS
1994  if (rField_is_Ring(r))
1995  {
1996    if ((ph!=NULL) && rField_has_Units(r))
1997    {
1998      number k = n_GetUnit(pGetCoeff(ph),r->cf);
1999      if (!n_IsOne(k,r->cf))
2000      {
2001        number tmpGMP = k;
2002        k = n_Invers(k,r->cf);
2003        n_Delete(&tmpGMP,r->cf);
2004        poly h = pNext(ph);
2005        p_SetCoeff(ph, n_Mult(pGetCoeff(ph), k,r->cf),r);
2006        while (h != NULL)
2007        {
2008          p_SetCoeff(h, n_Mult(pGetCoeff(h), k,r->cf),r);
2009          pIter(h);
2010        }
2011      }
2012      n_Delete(&k,r->cf);
2013    }
2014    return;
2015  }
2016#endif
2017  number h,d;
2018  poly p;
2019
2020  if(TEST_OPT_CONTENTSB) return;
2021  if(pNext(ph)==NULL)
2022  {
2023    p_SetCoeff(ph,n_Init(1,r->cf),r);
2024  }
2025  else
2026  {
2027    n_Normalize(pGetCoeff(ph),r->cf);
2028    if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2029    if (rField_is_Q(r))
2030    {
2031      h=p_InitContent(ph,r);
2032      p=ph;
2033    }
2034    else if (rField_is_Extension(r)
2035             &&
2036             (
2037              (rPar(r)>1) || rMinpolyIsNULL(r)
2038             )
2039            )
2040    {
2041      h=p_InitContent_a(ph,r);
2042      p=ph;
2043    }
2044    else
2045    {
2046      h=n_Copy(pGetCoeff(ph),r->cf);
2047      p = pNext(ph);
2048    }
2049    while (p!=NULL)
2050    {
2051      n_Normalize(pGetCoeff(p),r->cf);
2052      d=n_Gcd(h,pGetCoeff(p),r->cf);
2053      n_Delete(&h,r->cf);
2054      h = d;
2055      if(n_IsOne(h,r->cf))
2056      {
2057        break;
2058      }
2059      pIter(p);
2060    }
2061    p = ph;
2062    //number tmp;
2063    if(!n_IsOne(h,r->cf))
2064    {
2065      while (p!=NULL)
2066      {
2067        //d = nDiv(pGetCoeff(p),h);
2068        //tmp = nIntDiv(pGetCoeff(p),h);
2069        //if (!nEqual(d,tmp))
2070        //{
2071        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
2072        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
2073        //  nWrite(tmp);Print(StringAppendS("\n"));
2074        //}
2075        //nDelete(&tmp);
2076        d = n_IntDiv(pGetCoeff(p),h,r->cf);
2077        p_SetCoeff(p,d,r);
2078        pIter(p);
2079      }
2080    }
2081    n_Delete(&h,r->cf);
2082#ifdef HAVE_FACTORY
2083    if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
2084    {
2085      singclap_divide_content(ph, r);
2086      if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2087    }
2088#endif
2089    if (rField_is_Q_a(r))
2090    {
2091      //number hzz = nlInit(1, r->cf);
2092      h = nlInit(1, r->cf);
2093      p=ph;
2094      Werror("longalg missing");
2095      #if 0
2096      while (p!=NULL)
2097      { // each monom: coeff in Q_a
2098        lnumber c_n_n=(lnumber)pGetCoeff(p);
2099        poly c_n=c_n_n->z;
2100        while (c_n!=NULL)
2101        { // each monom: coeff in Q
2102          d=nlLcm(hzz,pGetCoeff(c_n),r->algring->cf);
2103          n_Delete(&hzz,r->algring->cf);
2104          hzz=d;
2105          pIter(c_n);
2106        }
2107        c_n=c_n_n->n;
2108        while (c_n!=NULL)
2109        { // each monom: coeff in Q
2110          d=nlLcm(h,pGetCoeff(c_n),r->algring->cf);
2111          n_Delete(&h,r->algring->cf);
2112          h=d;
2113          pIter(c_n);
2114        }
2115        pIter(p);
2116      }
2117      /* hzz contains the 1/lcm of all denominators in c_n_n->z*/
2118      /* h contains the 1/lcm of all denominators in c_n_n->n*/
2119      number htmp=nlInvers(h,r->algring->cf);
2120      number hzztmp=nlInvers(hzz,r->algring->cf);
2121      number hh=nlMult(hzz,h,r->algring->cf);
2122      nlDelete(&hzz,r->algring->cf);
2123      nlDelete(&h,r->algring->cf);
2124      number hg=nlGcd(hzztmp,htmp,r->algring->cf);
2125      nlDelete(&hzztmp,r->algring->cf);
2126      nlDelete(&htmp,r->algring->cf);
2127      h=nlMult(hh,hg,r->algring->cf);
2128      nlDelete(&hg,r->algring->cf);
2129      nlDelete(&hh,r->algring->cf);
2130      nlNormalize(h,r->algring->cf);
2131      if(!nlIsOne(h,r->algring->cf))
2132      {
2133        p=ph;
2134        while (p!=NULL)
2135        { // each monom: coeff in Q_a
2136          lnumber c_n_n=(lnumber)pGetCoeff(p);
2137          poly c_n=c_n_n->z;
2138          while (c_n!=NULL)
2139          { // each monom: coeff in Q
2140            d=nlMult(h,pGetCoeff(c_n),r->algring->cf);
2141            nlNormalize(d,r->algring->cf);
2142            nlDelete(&pGetCoeff(c_n),r->algring->cf);
2143            pGetCoeff(c_n)=d;
2144            pIter(c_n);
2145          }
2146          c_n=c_n_n->n;
2147          while (c_n!=NULL)
2148          { // each monom: coeff in Q
2149            d=nlMult(h,pGetCoeff(c_n),r->algring->cf);
2150            nlNormalize(d,r->algring->cf);
2151            nlDelete(&pGetCoeff(c_n),r->algring->cf);
2152            pGetCoeff(c_n)=d;
2153            pIter(c_n);
2154          }
2155          pIter(p);
2156        }
2157      }
2158      nlDelete(&h,r->algring->cf);
2159      #endif
2160    }
2161  }
2162}
2163#if 0 // currently not used
2164void p_SimpleContent(poly ph,int smax, const ring r)
2165{
2166  if(TEST_OPT_CONTENTSB) return;
2167  if (ph==NULL) return;
2168  if (pNext(ph)==NULL)
2169  {
2170    p_SetCoeff(ph,n_Init(1,r_cf),r);
2171    return;
2172  }
2173  if ((pNext(pNext(ph))==NULL)||(!rField_is_Q(r)))
2174  {
2175    return;
2176  }
2177  number d=p_InitContent(ph,r);
2178  if (nlSize(d,r->cf)<=smax)
2179  {
2180    //if (TEST_OPT_PROT) PrintS("G");
2181    return;
2182  }
2183  poly p=ph;
2184  number h=d;
2185  if (smax==1) smax=2;
2186  while (p!=NULL)
2187  {
2188#if 0
2189    d=nlGcd(h,pGetCoeff(p),r->cf);
2190    nlDelete(&h,r->cf);
2191    h = d;
2192#else
2193    nlInpGcd(h,pGetCoeff(p),r->cf);
2194#endif
2195    if(nlSize(h,r->cf)<smax)
2196    {
2197      //if (TEST_OPT_PROT) PrintS("g");
2198      return;
2199    }
2200    pIter(p);
2201  }
2202  p = ph;
2203  if (!nlGreaterZero(pGetCoeff(p),r->cf)) h=nlNeg(h,r->cf);
2204  if(nlIsOne(h,r->cf)) return;
2205  //if (TEST_OPT_PROT) PrintS("c");
2206  while (p!=NULL)
2207  {
2208#if 1
2209    d = nlIntDiv(pGetCoeff(p),h,r->cf);
2210    p_SetCoeff(p,d,r);
2211#else
2212    nlInpIntDiv(pGetCoeff(p),h,r->cf);
2213#endif
2214    pIter(p);
2215  }
2216  nlDelete(&h,r->cf);
2217}
2218#endif
2219
2220static number p_InitContent(poly ph, const ring r)
2221// only for coefficients in Q
2222#if 0
2223{
2224  assume(!TEST_OPT_CONTENTSB);
2225  assume(ph!=NULL);
2226  assume(pNext(ph)!=NULL);
2227  assume(rField_is_Q(r));
2228  if (pNext(pNext(ph))==NULL)
2229  {
2230    return nlGetNom(pGetCoeff(pNext(ph)),r->cf);
2231  }
2232  poly p=ph;
2233  number n1=nlGetNom(pGetCoeff(p),r->cf);
2234  pIter(p);
2235  number n2=nlGetNom(pGetCoeff(p),r->cf);
2236  pIter(p);
2237  number d;
2238  number t;
2239  loop
2240  {
2241    nlNormalize(pGetCoeff(p),r->cf);
2242    t=nlGetNom(pGetCoeff(p),r->cf);
2243    if (nlGreaterZero(t,r->cf))
2244      d=nlAdd(n1,t,r->cf);
2245    else
2246      d=nlSub(n1,t,r->cf);
2247    nlDelete(&t,r->cf);
2248    nlDelete(&n1,r->cf);
2249    n1=d;
2250    pIter(p);
2251    if (p==NULL) break;
2252    nlNormalize(pGetCoeff(p),r->cf);
2253    t=nlGetNom(pGetCoeff(p),r->cf);
2254    if (nlGreaterZero(t,r->cf))
2255      d=nlAdd(n2,t,r->cf);
2256    else
2257      d=nlSub(n2,t,r->cf);
2258    nlDelete(&t,r->cf);
2259    nlDelete(&n2,r->cf);
2260    n2=d;
2261    pIter(p);
2262    if (p==NULL) break;
2263  }
2264  d=nlGcd(n1,n2,r->cf);
2265  nlDelete(&n1,r->cf);
2266  nlDelete(&n2,r->cf);
2267  return d;
2268}
2269#else
2270{
2271  number d=pGetCoeff(ph);
2272  if(SR_HDL(d)&SR_INT) return d;
2273  int s=mpz_size1(d->z);
2274  int s2=-1;
2275  number d2;
2276  loop
2277  {
2278    pIter(ph);
2279    if(ph==NULL)
2280    {
2281      if (s2==-1) return nlCopy(d,r->cf);
2282      break;
2283    }
2284    if (SR_HDL(pGetCoeff(ph))&SR_INT)
2285    {
2286      s2=s;
2287      d2=d;
2288      s=0;
2289      d=pGetCoeff(ph);
2290      if (s2==0) break;
2291    }
2292    else
2293    if (mpz_size1((pGetCoeff(ph)->z))<=s)
2294    {
2295      s2=s;
2296      d2=d;
2297      d=pGetCoeff(ph);
2298      s=mpz_size1(d->z);
2299    }
2300  }
2301  return nlGcd(d,d2,r->cf);
2302}
2303#endif
2304
2305number p_InitContent_a(poly ph, const ring r)
2306// only for coefficients in K(a) anf K(a,...)
2307{
2308  number d=pGetCoeff(ph);
2309  int s=n_ParDeg(d,r->cf);
2310  if (s /* n_ParDeg(d)*/ <=1) return n_Copy(d,r->cf);
2311  int s2=-1;
2312  number d2;
2313  int ss;
2314  loop
2315  {
2316    pIter(ph);
2317    if(ph==NULL)
2318    {
2319      if (s2==-1) return n_Copy(d,r->cf);
2320      break;
2321    }
2322    if ((ss=n_ParDeg(pGetCoeff(ph),r->cf))<s)
2323    {
2324      s2=s;
2325      d2=d;
2326      s=ss;
2327      d=pGetCoeff(ph);
2328      if (s2<=1) break;
2329    }
2330  }
2331  return n_Gcd(d,d2,r->cf);
2332}
2333
2334
2335//void pContent(poly ph)
2336//{
2337//  number h,d;
2338//  poly p;
2339//
2340//  p = ph;
2341//  if(pNext(p)==NULL)
2342//  {
2343//    pSetCoeff(p,nInit(1));
2344//  }
2345//  else
2346//  {
2347//#ifdef PDEBUG
2348//    if (!pTest(p)) return;
2349//#endif
2350//    nNormalize(pGetCoeff(p));
2351//    if(!nGreaterZero(pGetCoeff(ph)))
2352//    {
2353//      ph = pNeg(ph);
2354//      nNormalize(pGetCoeff(p));
2355//    }
2356//    h=pGetCoeff(p);
2357//    pIter(p);
2358//    while (p!=NULL)
2359//    {
2360//      nNormalize(pGetCoeff(p));
2361//      if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p);
2362//      pIter(p);
2363//    }
2364//    h=nCopy(h);
2365//    p=ph;
2366//    while (p!=NULL)
2367//    {
2368//      d=n_Gcd(h,pGetCoeff(p));
2369//      nDelete(&h);
2370//      h = d;
2371//      if(nIsOne(h))
2372//      {
2373//        break;
2374//      }
2375//      pIter(p);
2376//    }
2377//    p = ph;
2378//    //number tmp;
2379//    if(!nIsOne(h))
2380//    {
2381//      while (p!=NULL)
2382//      {
2383//        d = nIntDiv(pGetCoeff(p),h);
2384//        pSetCoeff(p,d);
2385//        pIter(p);
2386//      }
2387//    }
2388//    nDelete(&h);
2389//#ifdef HAVE_FACTORY
2390//    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
2391//    {
2392//      pTest(ph);
2393//      singclap_divide_content(ph);
2394//      pTest(ph);
2395//    }
2396//#endif
2397//  }
2398//}
2399#if 0
2400void p_Content(poly ph, const ring r)
2401{
2402  number h,d;
2403  poly p;
2404
2405  if(pNext(ph)==NULL)
2406  {
2407    p_SetCoeff(ph,n_Init(1,r->cf),r);
2408  }
2409  else
2410  {
2411    n_Normalize(pGetCoeff(ph),r->cf);
2412    if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2413    h=n_Copy(pGetCoeff(ph),r->cf);
2414    p = pNext(ph);
2415    while (p!=NULL)
2416    {
2417      n_Normalize(pGetCoeff(p),r->cf);
2418      d=n_Gcd(h,pGetCoeff(p),r->cf);
2419      n_Delete(&h,r->cf);
2420      h = d;
2421      if(n_IsOne(h,r->cf))
2422      {
2423        break;
2424      }
2425      pIter(p);
2426    }
2427    p = ph;
2428    //number tmp;
2429    if(!n_IsOne(h,r->cf))
2430    {
2431      while (p!=NULL)
2432      {
2433        //d = nDiv(pGetCoeff(p),h);
2434        //tmp = nIntDiv(pGetCoeff(p),h);
2435        //if (!nEqual(d,tmp))
2436        //{
2437        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
2438        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
2439        //  nWrite(tmp);Print(StringAppendS("\n"));
2440        //}
2441        //nDelete(&tmp);
2442        d = n_IntDiv(pGetCoeff(p),h,r->cf);
2443        p_SetCoeff(p,d,r->cf);
2444        pIter(p);
2445      }
2446    }
2447    n_Delete(&h,r->cf);
2448#ifdef HAVE_FACTORY
2449    //if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
2450    //{
2451    //  singclap_divide_content(ph);
2452    //  if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r);
2453    //}
2454#endif
2455  }
2456}
2457#endif
2458/* ---------------------------------------------------------------------------*/
2459/* cleardenom suff                                                           */
2460poly p_Cleardenom(poly ph, const ring r)
2461{
2462  poly start=ph;
2463  number d, h;
2464  poly p;
2465
2466#ifdef HAVE_RINGS
2467  if (rField_is_Ring(r))
2468  {
2469    p_Content(ph,r);
2470    return start;
2471  }
2472#endif
2473  if (rField_is_Zp(r) && TEST_OPT_INTSTRATEGY) return start;
2474  p = ph;
2475  if(pNext(p)==NULL)
2476  {
2477    if (TEST_OPT_CONTENTSB)
2478    {
2479      number n=n_GetDenom(pGetCoeff(p),r->cf);
2480      if (!n_IsOne(n,r->cf))
2481      {
2482        number nn=n_Mult(pGetCoeff(p),n,r->cf);
2483        n_Normalize(nn,r->cf);
2484        p_SetCoeff(p,nn,r);
2485      }
2486      n_Delete(&n,r->cf);
2487    }
2488    else
2489      p_SetCoeff(p,n_Init(1,r->cf),r);
2490  }
2491  else
2492  {
2493    h = n_Init(1,r->cf);
2494    while (p!=NULL)
2495    {
2496      n_Normalize(pGetCoeff(p),r->cf);
2497      d=n_Lcm(h,pGetCoeff(p),r->cf);
2498      n_Delete(&h,r->cf);
2499      h=d;
2500      pIter(p);
2501    }
2502    /* contains the 1/lcm of all denominators */
2503    if(!n_IsOne(h,r->cf))
2504    {
2505      p = ph;
2506      while (p!=NULL)
2507      {
2508        /* should be:
2509        * number hh;
2510        * nGetDenom(p->coef,&hh);
2511        * nMult(&h,&hh,&d);
2512        * nNormalize(d);
2513        * nDelete(&hh);
2514        * nMult(d,p->coef,&hh);
2515        * nDelete(&d);
2516        * nDelete(&(p->coef));
2517        * p->coef =hh;
2518        */
2519        d=n_Mult(h,pGetCoeff(p),r->cf);
2520        n_Normalize(d,r->cf);
2521        p_SetCoeff(p,d,r);
2522        pIter(p);
2523      }
2524      n_Delete(&h,r->cf);
2525      if (n_GetChar(r->cf)==1)
2526      {
2527        loop
2528        {
2529          h = n_Init(1,r->cf);
2530          p=ph;
2531          while (p!=NULL)
2532          {
2533            d=n_Lcm(h,pGetCoeff(p),r->cf);
2534            n_Delete(&h,r->cf);
2535            h=d;
2536            pIter(p);
2537          }
2538          /* contains the 1/lcm of all denominators */
2539          if(!n_IsOne(h,r->cf))
2540          {
2541            p = ph;
2542            while (p!=NULL)
2543            {
2544              /* should be:
2545              * number hh;
2546              * nGetDenom(p->coef,&hh);
2547              * nMult(&h,&hh,&d);
2548              * nNormalize(d);
2549              * nDelete(&hh);
2550              * nMult(d,p->coef,&hh);
2551              * nDelete(&d);
2552              * nDelete(&(p->coef));
2553              * p->coef =hh;
2554              */
2555              d=n_Mult(h,pGetCoeff(p),r->cf);
2556              n_Normalize(d,r->cf);
2557              p_SetCoeff(p,d,r);
2558              pIter(p);
2559            }
2560            n_Delete(&h,r->cf);
2561          }
2562          else
2563          {
2564            n_Delete(&h,r->cf);
2565            break;
2566          }
2567        }
2568      }
2569    }
2570    if (h!=NULL) n_Delete(&h,r->cf);
2571
2572    p_Content(ph,r);
2573#ifdef HAVE_RATGRING
2574    if (rIsRatGRing(r))
2575    {
2576      /* quick unit detection in the rational case is done in gr_nc_bba */
2577      pContentRat(ph);
2578      start=ph;
2579    }
2580#endif
2581  }
2582  return start;
2583}
2584
2585void p_Cleardenom_n(poly ph,const ring r,number &c)
2586{
2587  number d, h;
2588  poly p;
2589
2590  p = ph;
2591  if(pNext(p)==NULL)
2592  {
2593    c=n_Invers(pGetCoeff(p),r->cf);
2594    p_SetCoeff(p,n_Init(1,r->cf),r);
2595  }
2596  else
2597  {
2598    h = n_Init(1,r->cf);
2599    while (p!=NULL)
2600    {
2601      n_Normalize(pGetCoeff(p),r->cf);
2602      d=n_Lcm(h,pGetCoeff(p),r->cf);
2603      n_Delete(&h,r->cf);
2604      h=d;
2605      pIter(p);
2606    }
2607    c=h;
2608    /* contains the 1/lcm of all denominators */
2609    if(!n_IsOne(h,r->cf))
2610    {
2611      p = ph;
2612      while (p!=NULL)
2613      {
2614        /* should be:
2615        * number hh;
2616        * nGetDenom(p->coef,&hh);
2617        * nMult(&h,&hh,&d);
2618        * nNormalize(d);
2619        * nDelete(&hh);
2620        * nMult(d,p->coef,&hh);
2621        * nDelete(&d);
2622        * nDelete(&(p->coef));
2623        * p->coef =hh;
2624        */
2625        d=n_Mult(h,pGetCoeff(p),r->cf);
2626        n_Normalize(d,r->cf);
2627        p_SetCoeff(p,d,r);
2628        pIter(p);
2629      }
2630      if (rField_is_Q_a(r))
2631      {
2632        loop
2633        {
2634          h = n_Init(1,r->cf);
2635          p=ph;
2636          while (p!=NULL)
2637          {
2638            d=n_Lcm(h,pGetCoeff(p),r->cf);
2639            n_Delete(&h,r->cf);
2640            h=d;
2641            pIter(p);
2642          }
2643          /* contains the 1/lcm of all denominators */
2644          if(!n_IsOne(h,r->cf))
2645          {
2646            p = ph;
2647            while (p!=NULL)
2648            {
2649              /* should be:
2650              * number hh;
2651              * nGetDenom(p->coef,&hh);
2652              * nMult(&h,&hh,&d);
2653              * nNormalize(d);
2654              * nDelete(&hh);
2655              * nMult(d,p->coef,&hh);
2656              * nDelete(&d);
2657              * nDelete(&(p->coef));
2658              * p->coef =hh;
2659              */
2660              d=n_Mult(h,pGetCoeff(p),r->cf);
2661              n_Normalize(d,r->cf);
2662              p_SetCoeff(p,d,r);
2663              pIter(p);
2664            }
2665            number t=n_Mult(c,h,r->cf);
2666            n_Delete(&c,r->cf);
2667            c=t;
2668          }
2669          else
2670          {
2671            break;
2672          }
2673          n_Delete(&h,r->cf);
2674        }
2675      }
2676    }
2677  }
2678}
2679
2680number p_GetAllDenom(poly ph, const ring r)
2681{
2682  number d=n_Init(1,r->cf);
2683  poly p = ph;
2684
2685  while (p!=NULL)
2686  {
2687    number h=n_GetDenom(pGetCoeff(p),r->cf);
2688    if (!n_IsOne(h,r->cf))
2689    {
2690      number dd=n_Mult(d,h,r->cf);
2691      n_Delete(&d,r->cf);
2692      d=dd;
2693    }
2694    n_Delete(&h,r->cf);
2695    pIter(p);
2696  }
2697  return d;
2698}
2699
2700int p_Size(poly p, const ring r)
2701{
2702  int count = 0;
2703  while ( p != NULL )
2704  {
2705    count+= n_Size( pGetCoeff( p ), r->cf );
2706    pIter( p );
2707  }
2708  return count;
2709}
2710
2711/*2
2712*make p homogeneous by multiplying the monomials by powers of x_varnum
2713*assume: deg(var(varnum))==1
2714*/
2715poly p_Homogen (poly p, int varnum, const ring r)
2716{
2717  pFDegProc deg;
2718  if (r->pLexOrder && (r->order[0]==ringorder_lp))
2719    deg=p_Totaldegree;
2720  else
2721    deg=r->pFDeg;
2722
2723  poly q=NULL, qn;
2724  int  o,ii;
2725  sBucket_pt bp;
2726
2727  if (p!=NULL)
2728  {
2729    if ((varnum < 1) || (varnum > rVar(r)))
2730    {
2731      return NULL;
2732    }
2733    o=deg(p,r);
2734    q=pNext(p);
2735    while (q != NULL)
2736    {
2737      ii=deg(q,r);
2738      if (ii>o) o=ii;
2739      pIter(q);
2740    }
2741    q = p_Copy(p,r);
2742    bp = sBucketCreate(r);
2743    while (q != NULL)
2744    {
2745      ii = o-deg(q,r);
2746      if (ii!=0)
2747      {
2748        p_AddExp(q,varnum, (long)ii,r);
2749        p_Setm(q,r);
2750      }
2751      qn = pNext(q);
2752      pNext(q) = NULL;
2753      sBucket_Add_p(bp, q, 1);
2754      q = qn;
2755    }
2756    sBucketDestroyAdd(bp, &q, &ii);
2757  }
2758  return q;
2759}
2760
2761/*2
2762*tests if p is homogeneous with respect to the actual weigths
2763*/
2764BOOLEAN p_IsHomogeneous (poly p, const ring r)
2765{
2766  poly qp=p;
2767  int o;
2768
2769  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
2770  pFDegProc d;
2771  if (r->pLexOrder && (r->order[0]==ringorder_lp))
2772    d=p_Totaldegree;
2773  else
2774    d=r->pFDeg;
2775  o = d(p,r);
2776  do
2777  {
2778    if (d(qp,r) != o) return FALSE;
2779    pIter(qp);
2780  }
2781  while (qp != NULL);
2782  return TRUE;
2783}
2784
2785/*----------utilities for syzygies--------------*/
2786BOOLEAN   p_VectorHasUnitB(poly p, int * k, const ring r)
2787{
2788  poly q=p,qq;
2789  int i;
2790
2791  while (q!=NULL)
2792  {
2793    if (p_LmIsConstantComp(q,r))
2794    {
2795      i = p_GetComp(q,r);
2796      qq = p;
2797      while ((qq != q) && (p_GetComp(qq,r) != i)) pIter(qq);
2798      if (qq == q)
2799      {
2800        *k = i;
2801        return TRUE;
2802      }
2803      else
2804        pIter(q);
2805    }
2806    else pIter(q);
2807  }
2808  return FALSE;
2809}
2810
2811void   p_VectorHasUnit(poly p, int * k, int * len, const ring r)
2812{
2813  poly q=p,qq;
2814  int i,j=0;
2815
2816  *len = 0;
2817  while (q!=NULL)
2818  {
2819    if (p_LmIsConstantComp(q,r))
2820    {
2821      i = p_GetComp(q,r);
2822      qq = p;
2823      while ((qq != q) && (p_GetComp(qq,r) != i)) pIter(qq);
2824      if (qq == q)
2825      {
2826       j = 0;
2827       while (qq!=NULL)
2828       {
2829         if (p_GetComp(qq,r)==i) j++;
2830        pIter(qq);
2831       }
2832       if ((*len == 0) || (j<*len))
2833       {
2834         *len = j;
2835         *k = i;
2836       }
2837      }
2838    }
2839    pIter(q);
2840  }
2841}
2842
2843poly p_TakeOutComp1(poly * p, int k, const ring r)
2844{
2845  poly q = *p;
2846
2847  if (q==NULL) return NULL;
2848
2849  poly qq=NULL,result = NULL;
2850
2851  if (p_GetComp(q,r)==k)
2852  {
2853    result = q; /* *p */
2854    while ((q!=NULL) && (p_GetComp(q,r)==k))
2855    {
2856      p_SetComp(q,0,r);
2857      p_SetmComp(q,r);
2858      qq = q;
2859      pIter(q);
2860    }
2861    *p = q;
2862    pNext(qq) = NULL;
2863  }
2864  if (q==NULL) return result;
2865//  if (pGetComp(q) > k) pGetComp(q)--;
2866  while (pNext(q)!=NULL)
2867  {
2868    if (p_GetComp(pNext(q),r)==k)
2869    {
2870      if (result==NULL)
2871      {
2872        result = pNext(q);
2873        qq = result;
2874      }
2875      else
2876      {
2877        pNext(qq) = pNext(q);
2878        pIter(qq);
2879      }
2880      pNext(q) = pNext(pNext(q));
2881      pNext(qq) =NULL;
2882      p_SetComp(qq,0,r);
2883      p_SetmComp(qq,r);
2884    }
2885    else
2886    {
2887      pIter(q);
2888//      if (pGetComp(q) > k) pGetComp(q)--;
2889    }
2890  }
2891  return result;
2892}
2893
2894poly p_TakeOutComp(poly * p, int k, const ring r)
2895{
2896  poly q = *p,qq=NULL,result = NULL;
2897
2898  if (q==NULL) return NULL;
2899  BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(r);
2900  if (p_GetComp(q,r)==k)
2901  {
2902    result = q;
2903    do
2904    {
2905      p_SetComp(q,0,r);
2906      if (use_setmcomp) p_SetmComp(q,r);
2907      qq = q;
2908      pIter(q);
2909    }
2910    while ((q!=NULL) && (p_GetComp(q,r)==k));
2911    *p = q;
2912    pNext(qq) = NULL;
2913  }
2914  if (q==NULL) return result;
2915  if (p_GetComp(q,r) > k)
2916  {
2917    p_SubComp(q,1,r);
2918    if (use_setmcomp) p_SetmComp(q,r);
2919  }
2920  poly pNext_q;
2921  while ((pNext_q=pNext(q))!=NULL)
2922  {
2923    if (p_GetComp(pNext_q,r)==k)
2924    {
2925      if (result==NULL)
2926      {
2927        result = pNext_q;
2928        qq = result;
2929      }
2930      else
2931      {
2932        pNext(qq) = pNext_q;
2933        pIter(qq);
2934      }
2935      pNext(q) = pNext(pNext_q);
2936      pNext(qq) =NULL;
2937      p_SetComp(qq,0,r);
2938      if (use_setmcomp) p_SetmComp(qq,r);
2939    }
2940    else
2941    {
2942      /*pIter(q);*/ q=pNext_q;
2943      if (p_GetComp(q,r) > k)
2944      {
2945        p_SubComp(q,1,r);
2946        if (use_setmcomp) p_SetmComp(q,r);
2947      }
2948    }
2949  }
2950  return result;
2951}
2952
2953// Splits *p into two polys: *q which consists of all monoms with
2954// component == comp and *p of all other monoms *lq == pLength(*q)
2955void p_TakeOutComp(poly *r_p, long comp, poly *r_q, int *lq, const ring r)
2956{
2957  spolyrec pp, qq;
2958  poly p, q, p_prev;
2959  int l = 0;
2960
2961#ifdef HAVE_ASSUME
2962  int lp = pLength(*r_p);
2963#endif
2964
2965  pNext(&pp) = *r_p;
2966  p = *r_p;
2967  p_prev = &pp;
2968  q = &qq;
2969
2970  while(p != NULL)
2971  {
2972    while (p_GetComp(p,r) == comp)
2973    {
2974      pNext(q) = p;
2975      pIter(q);
2976      p_SetComp(p, 0,r);
2977      p_SetmComp(p,r);
2978      pIter(p);
2979      l++;
2980      if (p == NULL)
2981      {
2982        pNext(p_prev) = NULL;
2983        goto Finish;
2984      }
2985    }
2986    pNext(p_prev) = p;
2987    p_prev = p;
2988    pIter(p);
2989  }
2990
2991  Finish:
2992  pNext(q) = NULL;
2993  *r_p = pNext(&pp);
2994  *r_q = pNext(&qq);
2995  *lq = l;
2996#ifdef HAVE_ASSUME
2997  assume(pLength(*r_p) + pLength(*r_q) == lp);
2998#endif
2999  p_Test(*r_p,r);
3000  p_Test(*r_q,r);
3001}
3002
3003void p_DeleteComp(poly * p,int k, const ring r)
3004{
3005  poly q;
3006
3007  while ((*p!=NULL) && (p_GetComp(*p,r)==k)) p_LmDelete(p,r);
3008  if (*p==NULL) return;
3009  q = *p;
3010  if (p_GetComp(q,r)>k)
3011  {
3012    p_SubComp(q,1,r);
3013    p_SetmComp(q,r);
3014  }
3015  while (pNext(q)!=NULL)
3016  {
3017    if (p_GetComp(pNext(q),r)==k)
3018      p_LmDelete(&(pNext(q)),r);
3019    else
3020    {
3021      pIter(q);
3022      if (p_GetComp(q,r)>k)
3023      {
3024        p_SubComp(q,1,r);
3025        p_SetmComp(q,r);
3026      }
3027    }
3028  }
3029}
3030
3031/*2
3032* convert a vector to a set of polys,
3033* allocates the polyset, (entries 0..(*len)-1)
3034* the vector will not be changed
3035*/
3036void  p_Vec2Polys(poly v, poly* *p, int *len, const ring r)
3037{
3038  poly h;
3039  int k;
3040
3041  *len=p_MaxComp(v,r);
3042  if (*len==0) *len=1;
3043  *p=(poly*)omAlloc0((*len)*sizeof(poly));
3044  while (v!=NULL)
3045  {
3046    h=p_Head(v,r);
3047    k=p_GetComp(h,r);
3048    p_SetComp(h,0,r);
3049    (*p)[k-1]=p_Add_q((*p)[k-1],h,r);
3050    pIter(v);
3051  }
3052}
3053
3054/* -------------------------------------------------------- */
3055/*2
3056* change all global variables to fit the description of the new ring
3057*/
3058
3059void p_SetGlobals(const ring r, BOOLEAN complete)
3060{
3061  if (r->ppNoether!=NULL) p_Delete(&r->ppNoether,r);
3062
3063  if (complete)
3064  {
3065    test &= ~ TEST_RINGDEP_OPTS;
3066    test |= r->options;
3067  }
3068}
3069//
3070// resets the pFDeg and pLDeg: if pLDeg is not given, it is
3071// set to currRing->pLDegOrig, i.e. to the respective LDegProc which
3072// only uses pFDeg (and not p_Deg, or pTotalDegree, etc)
3073void pSetDegProcs(ring r, pFDegProc new_FDeg, pLDegProc new_lDeg)
3074{
3075  assume(new_FDeg != NULL);
3076  r->pFDeg = new_FDeg;
3077
3078  if (new_lDeg == NULL)
3079    new_lDeg = r->pLDegOrig;
3080
3081  r->pLDeg = new_lDeg;
3082}
3083
3084// restores pFDeg and pLDeg:
3085void pRestoreDegProcs(ring r, pFDegProc old_FDeg, pLDegProc old_lDeg)
3086{
3087  assume(old_FDeg != NULL && old_lDeg != NULL);
3088  r->pFDeg = old_FDeg;
3089  r->pLDeg = old_lDeg;
3090}
3091
3092/*-------- several access procedures to monomials -------------------- */
3093/*
3094* the module weights for std
3095*/
3096static pFDegProc pOldFDeg;
3097static pLDegProc pOldLDeg;
3098static intvec * pModW;
3099static BOOLEAN pOldLexOrder;
3100
3101static long pModDeg(poly p, ring r)
3102{
3103  long d=pOldFDeg(p, r);
3104  int c=p_GetComp(p, r);
3105  if ((c>0) && ((r->pModW)->range(c-1))) d+= (*(r->pModW))[c-1];
3106  return d;
3107  //return pOldFDeg(p, r)+(*pModW)[p_GetComp(p, r)-1];
3108}
3109
3110void p_SetModDeg(intvec *w, ring r)
3111{
3112  if (w!=NULL)
3113  {
3114    r->pModW = w;
3115    pOldFDeg = r->pFDeg;
3116    pOldLDeg = r->pLDeg;
3117    pOldLexOrder = r->pLexOrder;
3118    pSetDegProcs(r,pModDeg);
3119    r->pLexOrder = TRUE;
3120  }
3121  else
3122  {
3123    r->pModW = NULL;
3124    pRestoreDegProcs(r,pOldFDeg, pOldLDeg);
3125    r->pLexOrder = pOldLexOrder;
3126  }
3127}
3128
3129/*2
3130* handle memory request for sets of polynomials (ideals)
3131* l is the length of *p, increment is the difference (may be negative)
3132*/
3133void pEnlargeSet(poly* *p, int l, int increment)
3134{
3135  poly* h;
3136
3137  h=(poly*)omReallocSize((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly));
3138  if (increment>0)
3139  {
3140    //for (i=l; i<l+increment; i++)
3141    //  h[i]=NULL;
3142    memset(&(h[l]),0,increment*sizeof(poly));
3143  }
3144  *p=h;
3145}
3146
3147/*2
3148*divides p1 by its leading coefficient
3149*/
3150void p_Norm(poly p1, const ring r)
3151{
3152#ifdef HAVE_RINGS
3153  if (rField_is_Ring(r))
3154  {
3155    if (!n_IsUnit(pGetCoeff(p1), r->cf)) return;
3156    // Werror("p_Norm not possible in the case of coefficient rings.");
3157  }
3158  else
3159#endif
3160  if (p1!=NULL)
3161  {
3162    if (pNext(p1)==NULL)
3163    {
3164      p_SetCoeff(p1,n_Init(1,r->cf),r);
3165      return;
3166    }
3167    poly h;
3168    if (!n_IsOne(pGetCoeff(p1),r->cf))
3169    {
3170      number k, c;
3171      n_Normalize(pGetCoeff(p1),r->cf);
3172      k = pGetCoeff(p1);
3173      c = n_Init(1,r->cf);
3174      pSetCoeff0(p1,c);
3175      h = pNext(p1);
3176      while (h!=NULL)
3177      {
3178        c=n_Div(pGetCoeff(h),k,r->cf);
3179        // no need to normalize: Z/p, R
3180        // normalize already in nDiv: Q_a, Z/p_a
3181        // remains: Q
3182        if (rField_is_Q(r) && (!n_IsOne(c,r->cf))) n_Normalize(c,r->cf);
3183        p_SetCoeff(h,c,r);
3184        pIter(h);
3185      }
3186      n_Delete(&k,r->cf);
3187    }
3188    else
3189    {
3190      //if (r->cf->cfNormalize != nDummy2) //TODO: OPTIMIZE
3191      {
3192        h = pNext(p1);
3193        while (h!=NULL)
3194        {
3195          n_Normalize(pGetCoeff(h),r->cf);
3196          pIter(h);
3197        }
3198      }
3199    }
3200  }
3201}
3202
3203/*2
3204*normalize all coefficients
3205*/
3206void p_Normalize(poly p,const ring r)
3207{
3208  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
3209  while (p!=NULL)
3210  {
3211#ifdef LDEBUG
3212    n_Test(pGetCoeff(p), r->cf);
3213#endif
3214    n_Normalize(pGetCoeff(p),r->cf);
3215    pIter(p);
3216  }
3217}
3218
3219// splits p into polys with Exp(n) == 0 and Exp(n) != 0
3220// Poly with Exp(n) != 0 is reversed
3221static void p_SplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero, const ring r)
3222{
3223  if (p == NULL)
3224  {
3225    *non_zero = NULL;
3226    *zero = NULL;
3227    return;
3228  }
3229  spolyrec sz;
3230  poly z, n_z, next;
3231  z = &sz;
3232  n_z = NULL;
3233
3234  while(p != NULL)
3235  {
3236    next = pNext(p);
3237    if (p_GetExp(p, n,r) == 0)
3238    {
3239      pNext(z) = p;
3240      pIter(z);
3241    }
3242    else
3243    {
3244      pNext(p) = n_z;
3245      n_z = p;
3246    }
3247    p = next;
3248  }
3249  pNext(z) = NULL;
3250  *zero = pNext(&sz);
3251  *non_zero = n_z;
3252}
3253/*3
3254* substitute the n-th variable by 1 in p
3255* destroy p
3256*/
3257static poly p_Subst1 (poly p,int n, const ring r)
3258{
3259  poly qq=NULL, result = NULL;
3260  poly zero=NULL, non_zero=NULL;
3261
3262  // reverse, so that add is likely to be linear
3263  p_SplitAndReversePoly(p, n, &non_zero, &zero,r);
3264
3265  while (non_zero != NULL)
3266  {
3267    assume(p_GetExp(non_zero, n,r) != 0);
3268    qq = non_zero;
3269    pIter(non_zero);
3270    qq->next = NULL;
3271    p_SetExp(qq,n,0,r);
3272    p_Setm(qq,r);
3273    result = p_Add_q(result,qq,r);
3274  }
3275  p = p_Add_q(result, zero,r);
3276  p_Test(p,r);
3277  return p;
3278}
3279
3280/*3
3281* substitute the n-th variable by number e in p
3282* destroy p
3283*/
3284static poly p_Subst2 (poly p,int n, number e, const ring r)
3285{
3286  assume( ! n_IsZero(e,r->cf) );
3287  poly qq,result = NULL;
3288  number nn, nm;
3289  poly zero, non_zero;
3290
3291  // reverse, so that add is likely to be linear
3292  p_SplitAndReversePoly(p, n, &non_zero, &zero,r);
3293
3294  while (non_zero != NULL)
3295  {
3296    assume(p_GetExp(non_zero, n, r) != 0);
3297    qq = non_zero;
3298    pIter(non_zero);
3299    qq->next = NULL;
3300    n_Power(e, p_GetExp(qq, n, r), &nn,r->cf);
3301    nm = n_Mult(nn, pGetCoeff(qq),r->cf);
3302#ifdef HAVE_RINGS
3303    if (n_IsZero(nm,r->cf))
3304    {
3305      p_LmFree(&qq,r);
3306      n_Delete(&nm,r->cf);
3307    }
3308    else
3309#endif
3310    {
3311      p_SetCoeff(qq, nm,r);
3312      p_SetExp(qq, n, 0,r);
3313      p_Setm(qq,r);
3314      result = p_Add_q(result,qq,r);
3315    }
3316    n_Delete(&nn,r->cf);
3317  }
3318  p = p_Add_q(result, zero,r);
3319  p_Test(p,r);
3320  return p;
3321}
3322
3323
3324/* delete monoms whose n-th exponent is different from zero */
3325static poly p_Subst0(poly p, int n, const ring r)
3326{
3327  spolyrec res;
3328  poly h = &res;
3329  pNext(h) = p;
3330
3331  while (pNext(h)!=NULL)
3332  {
3333    if (p_GetExp(pNext(h),n,r)!=0)
3334    {
3335      p_LmDelete(&pNext(h),r);
3336    }
3337    else
3338    {
3339      pIter(h);
3340    }
3341  }
3342  p_Test(pNext(&res),r);
3343  return pNext(&res);
3344}
3345
3346/*2
3347* substitute the n-th variable by e in p
3348* destroy p
3349*/
3350poly p_Subst(poly p, int n, poly e, const ring r)
3351{
3352  if (e == NULL) return p_Subst0(p, n,r);
3353
3354  if (p_IsConstant(e,r))
3355  {
3356    if (n_IsOne(pGetCoeff(e),r->cf)) return p_Subst1(p,n,r);
3357    else return p_Subst2(p, n, pGetCoeff(e),r);
3358  }
3359
3360#ifdef HAVE_PLURAL
3361  if (rIsPluralRing(r))
3362  {
3363    return nc_pSubst(p,n,e,r);
3364  }
3365#endif
3366
3367  int exponent,i;
3368  poly h, res, m;
3369  int *me,*ee;
3370  number nu,nu1;
3371
3372  me=(int *)omAlloc((rVar(r)+1)*sizeof(int));
3373  ee=(int *)omAlloc((rVar(r)+1)*sizeof(int));
3374  if (e!=NULL) p_GetExpV(e,ee,r);
3375  res=NULL;
3376  h=p;
3377  while (h!=NULL)
3378  {
3379    if ((e!=NULL) || (p_GetExp(h,n,r)==0))
3380    {
3381      m=p_Head(h,r);
3382      p_GetExpV(m,me,r);
3383      exponent=me[n];
3384      me[n]=0;
3385      for(i=rVar(r);i>0;i--)
3386        me[i]+=exponent*ee[i];
3387      p_SetExpV(m,me,r);
3388      if (e!=NULL)
3389      {
3390        n_Power(pGetCoeff(e),exponent,&nu,r->cf);
3391        nu1=n_Mult(pGetCoeff(m),nu,r->cf);
3392        n_Delete(&nu,r->cf);
3393        p_SetCoeff(m,nu1,r);
3394      }
3395      res=p_Add_q(res,m,r);
3396    }
3397    p_LmDelete(&h,r);
3398  }
3399  omFreeSize((ADDRESS)me,(rVar(r)+1)*sizeof(int));
3400  omFreeSize((ADDRESS)ee,(rVar(r)+1)*sizeof(int));
3401  return res;
3402}
3403/*2
3404*returns a re-ordered copy of a polynomial, with permutation of the variables
3405*/
3406poly p_PermPoly (poly p, int * perm, const ring oldRing, const ring dst,
3407       nMapFunc nMap, int *par_perm, int OldPar)
3408{
3409  int OldpVariables = oldRing->N;
3410  poly result = NULL;
3411  poly result_last = NULL;
3412  poly aq=NULL; /* the map coefficient */
3413  poly qq; /* the mapped monomial */
3414
3415  while (p != NULL)
3416  {
3417    if ((OldPar==0)||(rField_is_GF(oldRing)))
3418    {
3419      qq = p_Init(dst);
3420      number n=nMap(pGetCoeff(p),oldRing->cf,dst->cf);
3421      if ((!rMinpolyIsNULL(dst))
3422      && ((rField_is_Zp_a(dst)) || (rField_is_Q_a(dst))))
3423      {
3424        n_Normalize(n,dst->cf);
3425      }
3426      pGetCoeff(qq)=n;
3427    // coef may be zero:  pTest(qq);
3428    }
3429    else
3430    {
3431      qq=p_One(dst);
3432      WerrorS("longalg missing");
3433      #if 0
3434      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
3435      if ((!rMinpolyIsNULL(dst))
3436      && ((rField_is_Zp_a(dst)) || (rField_is_Q_a(dst))))
3437      {
3438        poly tmp=aq;
3439        while (tmp!=NULL)
3440        {
3441          number n=pGetCoeff(tmp);
3442          n_Normalize(n,dst->cf);
3443          pGetCoeff(tmp)=n;
3444          pIter(tmp);
3445        }
3446      }
3447      p_Test(aq,dst);
3448      #endif
3449    }
3450    if (rRing_has_Comp(dst)) p_SetComp(qq, p_GetComp(p,oldRing),dst);
3451    if (n_IsZero(pGetCoeff(qq),dst->cf))
3452    {
3453      p_LmDelete(&qq,dst);
3454    }
3455    else
3456    {
3457      int i;
3458      int mapped_to_par=0;
3459      for(i=1; i<=OldpVariables; i++)
3460      {
3461        int e=p_GetExp(p,i,oldRing);
3462        if (e!=0)
3463        {
3464          if (perm==NULL)
3465          {
3466            p_SetExp(qq,i, e, dst);
3467          }
3468          else if (perm[i]>0)
3469            p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, dst);
3470          else if (perm[i]<0)
3471          {
3472            if (rField_is_GF(dst))
3473            {
3474              number c=pGetCoeff(qq);
3475              number ee=n_Par(1,dst->cf);
3476              number eee;n_Power(ee,e,&eee,dst->cf); //nfDelete(ee,dst);
3477              ee=n_Mult(c,eee,dst->cf);
3478              //nfDelete(c,dst);nfDelete(eee,dst);
3479              pSetCoeff0(qq,ee);
3480            }
3481            else
3482            {
3483              WerrorS("longalg missing");
3484              #if 0
3485              lnumber c=(lnumber)pGetCoeff(qq);
3486              if (c->z->next==NULL)
3487                p_AddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/,dst->algring);
3488              else /* more difficult: we have really to multiply: */
3489              {
3490                lnumber mmc=(lnumber)naInit(1,dst);
3491                p_SetExp(mmc->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/,dst->algring);
3492                p_Setm(mmc->z,dst->algring->cf);
3493                pGetCoeff(qq)=n_Mult((number)c,(number)mmc,dst->cf);
3494                n_Delete((number *)&c,dst->cf);
3495                n_Delete((number *)&mmc,dst->cf);
3496              }
3497              mapped_to_par=1;
3498              #endif
3499            }
3500          }
3501          else
3502          {
3503            /* this variable maps to 0 !*/
3504            p_LmDelete(&qq,dst);
3505            break;
3506          }
3507        }
3508      }
3509      if (mapped_to_par
3510      && (!rMinpolyIsNULL(dst)))
3511      {
3512        number n=pGetCoeff(qq);
3513        n_Normalize(n,dst->cf);
3514        pGetCoeff(qq)=n;
3515      }
3516    }
3517    pIter(p);
3518#if 1
3519    if (qq!=NULL)
3520    {
3521      p_Setm(qq,dst);
3522      p_Test(aq,dst);
3523      p_Test(qq,dst);
3524      if (aq!=NULL) qq=p_Mult_q(aq,qq,dst);
3525      aq = qq;
3526      while (pNext(aq) != NULL) pIter(aq);
3527      if (result_last==NULL)
3528      {
3529        result=qq;
3530      }
3531      else
3532      {
3533        pNext(result_last)=qq;
3534      }
3535      result_last=aq;
3536      aq = NULL;
3537    }
3538    else if (aq!=NULL)
3539    {
3540      p_Delete(&aq,dst);
3541    }
3542  }
3543  result=p_SortAdd(result,dst);
3544#else
3545  //  if (qq!=NULL)
3546  //  {
3547  //    pSetm(qq);
3548  //    pTest(qq);
3549  //    pTest(aq);
3550  //    if (aq!=NULL) qq=pMult(aq,qq);
3551  //    aq = qq;
3552  //    while (pNext(aq) != NULL) pIter(aq);
3553  //    pNext(aq) = result;
3554  //    aq = NULL;
3555  //    result = qq;
3556  //  }
3557  //  else if (aq!=NULL)
3558  //  {
3559  //    pDelete(&aq);
3560  //  }
3561  //}
3562  //p = result;
3563  //result = NULL;
3564  //while (p != NULL)
3565  //{
3566  //  qq = p;
3567  //  pIter(p);
3568  //  qq->next = NULL;
3569  //  result = pAdd(result, qq);
3570  //}
3571#endif
3572  p_Test(result,dst);
3573  return result;
3574}
3575/**************************************************************
3576 *
3577 * Jet
3578 *
3579 **************************************************************/
3580
3581poly pp_Jet(poly p, int m, const ring R)
3582{
3583  poly r=NULL;
3584  poly t=NULL;
3585
3586  while (p!=NULL)
3587  {
3588    if (p_Totaldegree(p,R)<=m)
3589    {
3590      if (r==NULL)
3591        r=p_Head(p,R);
3592      else
3593      if (t==NULL)
3594      {
3595        pNext(r)=p_Head(p,R);
3596        t=pNext(r);
3597      }
3598      else
3599      {
3600        pNext(t)=p_Head(p,R);
3601        pIter(t);
3602      }
3603    }
3604    pIter(p);
3605  }
3606  return r;
3607}
3608
3609poly p_Jet(poly p, int m,const ring R)
3610{
3611  poly t=NULL;
3612
3613  while((p!=NULL) && (p_Totaldegree(p,R)>m)) p_LmDelete(&p,R);
3614  if (p==NULL) return NULL;
3615  poly r=p;
3616  while (pNext(p)!=NULL)
3617  {
3618    if (p_Totaldegree(pNext(p),R)>m)
3619    {
3620      p_LmDelete(&pNext(p),R);
3621    }
3622    else
3623      pIter(p);
3624  }
3625  return r;
3626}
3627
3628poly pp_JetW(poly p, int m, short *w, const ring R)
3629{
3630  poly r=NULL;
3631  poly t=NULL;
3632  while (p!=NULL)
3633  {
3634    if (totaldegreeWecart_IV(p,R,w)<=m)
3635    {
3636      if (r==NULL)
3637        r=p_Head(p,R);
3638      else
3639      if (t==NULL)
3640      {
3641        pNext(r)=p_Head(p,R);
3642        t=pNext(r);
3643      }
3644      else
3645      {
3646        pNext(t)=p_Head(p,R);
3647        pIter(t);
3648      }
3649    }
3650    pIter(p);
3651  }
3652  return r;
3653}
3654
3655poly p_JetW(poly p, int m, short *w, const ring R)
3656{
3657  while((p!=NULL) && (totaldegreeWecart_IV(p,R,w)>m)) p_LmDelete(&p,R);
3658  if (p==NULL) return NULL;
3659  poly r=p;
3660  while (pNext(p)!=NULL)
3661  {
3662    if (totaldegreeWecart_IV(pNext(p),R,w)>m)
3663    {
3664      p_LmDelete(&pNext(p),R);
3665    }
3666    else
3667      pIter(p);
3668  }
3669  return r;
3670}
3671
3672/*************************************************************/
3673int p_MinDeg(poly p,intvec *w, const ring R)
3674{
3675  if(p==NULL)
3676    return -1;
3677  int d=-1;
3678  while(p!=NULL)
3679  {
3680    int d0=0;
3681    for(int j=0;j<rVar(R);j++)
3682      if(w==NULL||j>=w->length())
3683        d0+=p_GetExp(p,j+1,R);
3684      else
3685        d0+=(*w)[j]*p_GetExp(p,j+1,R);
3686    if(d0<d||d==-1)
3687      d=d0;
3688    pIter(p);
3689  }
3690  return d;
3691}
3692
3693/***************************************************************/
3694
3695poly p_Series(int n,poly p,poly u, intvec *w, const ring R)
3696{
3697  short *ww=iv2array(w,R);
3698  if(p!=NULL)
3699  {
3700    if(u==NULL)
3701      p=p_JetW(p,n,ww,R);
3702    else
3703      p=p_JetW(p_Mult_q(p,p_Invers(n-p_MinDeg(p,w,R),u,w,R),R),n,ww,R);
3704  }
3705  omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short));
3706  return p;
3707}
3708
3709poly p_Invers(int n,poly u,intvec *w, const ring R)
3710{
3711  if(n<0)
3712    return NULL;
3713  number u0=n_Invers(pGetCoeff(u),R->cf);
3714  poly v=p_NSet(u0,R);
3715  if(n==0)
3716    return v;
3717  short *ww=iv2array(w,R);
3718  poly u1=p_JetW(p_Sub(p_One(R),p_Mult_nn(u,u0,R),R),n,ww,R);
3719  if(u1==NULL)
3720  {
3721    omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short));
3722    return v;
3723  }
3724  poly v1=p_Mult_nn(p_Copy(u1,R),u0,R);
3725  v=p_Add_q(v,p_Copy(v1,R),R);
3726  for(int i=n/p_MinDeg(u1,w,R);i>1;i--)
3727  {
3728    v1=p_JetW(p_Mult_q(v1,p_Copy(u1,R),R),n,ww,R);
3729    v=p_Add_q(v,p_Copy(v1,R),R);
3730  }
3731  p_Delete(&u1,R);
3732  p_Delete(&v1,R);
3733  omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short));
3734  return v;
3735}
3736
3737BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r)
3738{
3739  while ((p1 != NULL) && (p2 != NULL))
3740  {
3741    if (! p_LmEqual(p1, p2,r))
3742      return FALSE;
3743    if (! n_Equal(p_GetCoeff(p1,r), p_GetCoeff(p2,r),r ))
3744      return FALSE;
3745    pIter(p1);
3746    pIter(p2);
3747  }
3748  return (p1==p2);
3749}
3750
3751/*2
3752*returns TRUE if p1 is a skalar multiple of p2
3753*assume p1 != NULL and p2 != NULL
3754*/
3755BOOLEAN p_ComparePolys(poly p1,poly p2, const ring r)
3756{
3757  number n,nn;
3758  pAssume(p1 != NULL && p2 != NULL);
3759
3760  if (!p_LmEqual(p1,p2,r)) //compare leading mons
3761      return FALSE;
3762  if ((pNext(p1)==NULL) && (pNext(p2)!=NULL))
3763     return FALSE;
3764  if ((pNext(p2)==NULL) && (pNext(p1)!=NULL))
3765     return FALSE;
3766  if (pLength(p1) != pLength(p2))
3767    return FALSE;
3768#ifdef HAVE_RINGS
3769  if (rField_is_Ring(r))
3770  {
3771    if (!n_DivBy(p_GetCoeff(p1, r), p_GetCoeff(p2, r), r)) return FALSE;
3772  }
3773#endif
3774  n=n_Div(p_GetCoeff(p1,r),p_GetCoeff(p2,r),r);
3775  while ((p1 != NULL) /*&& (p2 != NULL)*/)
3776  {
3777    if ( ! p_LmEqual(p1, p2,r))
3778    {
3779        n_Delete(&n, r);
3780        return FALSE;
3781    }
3782    if (!n_Equal(p_GetCoeff(p1, r), nn = n_Mult(p_GetCoeff(p2, r),n, r), r))
3783    {
3784      n_Delete(&n, r);
3785      n_Delete(&nn, r);
3786      return FALSE;
3787    }
3788    n_Delete(&nn, r);
3789    pIter(p1);
3790    pIter(p2);
3791  }
3792  n_Delete(&n, r);
3793  return TRUE;
3794}
3795
3796
3797/***************************************************************
3798 *
3799 * p_ShallowDelete
3800 *
3801 ***************************************************************/
3802#undef LINKAGE
3803#define LINKAGE
3804#undef p_Delete__T
3805#define p_Delete__T p_ShallowDelete
3806#undef n_Delete__T
3807#define n_Delete__T(n, r) ((void)0)
3808
3809#include <polys/templates/p_Delete__T.cc>
3810
3811#ifdef HAVE_RINGS
3812/* TRUE iff LT(f) | LT(g) */
3813BOOLEAN p_DivisibleByRingCase(poly f, poly g, const ring r)
3814{
3815  int exponent;
3816  for(int i = (int)r->N; i; i--)
3817  {
3818    exponent = p_GetExp(g, i, r) - p_GetExp(f, i, r);
3819    if (exponent < 0) return FALSE;
3820  }
3821  return n_DivBy(p_GetCoeff(g,r), p_GetCoeff(f,r),r->cf);
3822}
3823#endif
Note: See TracBrowser for help on using the repository browser.