source: git/kernel/ratgring.cc @ 61944d0

spielwiese
Last change on this file since 61944d0 was 14db39, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: pConentRat git-svn-id: file:///usr/local/Singular/svn/trunk@11468 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 19.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    ratgring.cc
6 *  Purpose: Ore-noncommutative kernel procedures
7 *  Author:  levandov (Viktor Levandovsky)
8 *  Created: 8/00 - 11/00
9 *  Version: $Id: ratgring.cc,v 1.17 2009-02-25 17:01:00 Singular Exp $
10 *******************************************************************/
11#include "mod2.h"
12#include "ratgring.h"
13#ifdef HAVE_RATGRING
14#include "gring.h"
15#include "febase.h"
16#include "ring.h"
17#include "polys.h"
18#include "numbers.h"
19#include "ideals.h"
20#include "matpol.h"
21#include "kbuckets.h"
22#include "kstd1.h"
23#include "sbuckets.h"
24#include "prCopy.h"
25#include "p_Mult_q.h"
26#include "clapsing.h"
27
28void pLcmRat(poly a, poly b, poly m, int rat_shift)
29{
30  /* rat_shift is the last exp one should count with */
31  int i;
32  for (i=pVariables; i>=rat_shift; i--)
33  {
34    pSetExp(m,i, si_max( pGetExp(a,i), pGetExp(b,i)));
35  }
36  pSetComp(m, si_max(pGetComp(a), pGetComp(b)));
37  /* Don't do a pSetm here, otherwise hres/lres chockes */ 
38}
39
40/*2
41* returns the rational LCM of the head terms of a and b
42* without coefficient!!!
43*/
44poly p_LcmRat(const poly a, const poly b, const long lCompM, const ring r)
45{
46  poly m = // p_One( r);
47          p_Init(r);
48
49  const int pVariables = r->N;
50
51  //  for (int i = pVariables; i>=r->real_var_start; i--)
52  for (int i = r->real_var_end; i>=r->real_var_start; i--)
53  {
54    const int lExpA = p_GetExp (a, i, r);
55    const int lExpB = p_GetExp (b, i, r);
56
57    p_SetExp (m, i, si_max(lExpA, lExpB), r);
58  }
59
60  p_SetComp (m, lCompM, r);
61  p_Setm(m,r);
62  n_New(&(p_GetCoeff(m, r)), r);
63
64  return(m);
65};
66
67// void pLcmRat(poly a, poly b, poly m, poly pshift)
68// {
69//   /* shift is the exp of rational elements */
70//   int i;
71//   for (i=pVariables; i; i--)
72//   {
73//     if (!pGetExp(pshift,i))
74//     {
75//       pSetExp(m,i, si_max( pGetExp(a,i), pGetExp(b,i)));
76//     }
77//     else
78//     {
79//       /* do we really need it? */
80//       pSetExp(m,i,0);
81//     }
82//   }
83//   pSetComp(m, si_max(pGetComp(a), pGetComp(b)));
84//   /* Don't do a pSetm here, otherwise hres/lres chockes */
85// }
86
87/* returns a subpoly of p, s.t. its monomials have the same D-part */
88
89poly p_HeadRat(poly p, int ishift, ring r)
90{
91  poly q   = pNext(p);
92  if (q == NULL) return p;
93  poly res = p_Head(p,r);
94  while ( (q!=NULL) && (p_Comp_k_n(p, q, ishift+1, r)))
95  {
96    res = p_Add_q(res,p_Head(q,r),r);
97    q   = pNext(q);
98  }
99  return res;
100}
101
102/* TO TEST!!! */
103/* returns x-coeff of p, i.e. a poly in x, s.t. corresponding xd-monomials
104have the same D-part
105does not destroy p
106*/
107
108poly p_GetCoeffRat(poly p, int ishift, ring r)
109{
110  poly q   = pNext(p);
111  poly res; // = p_Head(p,r);
112  res = p_GetExp_k_n(p, ishift+1, r->N, r);
113  p_SetCoeff(res,n_Copy(p_GetCoeff(p,r),r),r);
114  poly s;
115  while ((q!= NULL) && (p_Comp_k_n(p, q, ishift+1, r)))
116  {
117    s   = p_GetExp_k_n(q, ishift+1, r->N, r);
118    p_SetCoeff(s,n_Copy(p_GetCoeff(q,r),r),r);
119    res = p_Add_q(res,s,r);
120    q   = pNext(q);
121  }
122  return res;
123}
124
125void p_LmDeleteAndNextRat(poly *p, int ishift, ring r)
126{
127  /* modifies p*/
128  //  Print("start: "); Print(" "); p_wrp(*p,r);
129  p_LmCheckPolyRing2(*p, r);
130  poly q = p_Head(*p,r);
131  // in the next line ishift is correct
132  while ( ( (*p)!=NULL ) && ( p_Comp_k_n(*p, q, ishift, r) ))
133  {
134    p_LmDelete(p,r);
135    //    Print("while: ");p_wrp(*p,r);Print(" ");
136  }
137  //  p_wrp(*p,r);Print(" ");
138  //  PrintS("end\n");
139  p_LmDelete(&q,r);
140}
141
142/* to test!!! */
143/* ExpVector(pr) = ExpVector(p1) - ExpVector(p2) */
144void p_ExpVectorDiffRat(poly pr, poly p1, poly p2, int ishift, ring r)
145{
146  p_LmCheckPolyRing1(p1, r);
147  p_LmCheckPolyRing1(p2, r);
148  p_LmCheckPolyRing1(pr, r); 
149  int i;
150  poly t=pr;
151  Exponent_t e1,e2;
152  for (i=ishift+1; i<=r->N; i++)
153  {
154    e1 = p_GetExp(p1, i, r);
155    e2 = p_GetExp(p2, i, r);
156    //    pAssume1(p_GetExp(p1, i, r) >= p_GetExp(p2, i, r));
157    if (e1 < e2)
158    {
159#ifdef PDEBUG
160      Print("negative ExpVectorDiff\n");
161#endif   
162      p_Delete(&t,r);
163      break;
164    }
165    else
166    {
167      p_SetExp(t,i, e1-e2,r);
168    }
169  }
170  p_Setm(t,r);
171}
172
173/* returns ideal (u,v) s.t. up + vq = 0 */
174
175ideal ncGCD2(poly p, poly q, const ring r)
176{
177  // todo: must destroy p,q
178  intvec *w = NULL;
179  ideal h = idInit(2,1);
180  h->m[0] = p_Copy(p,r);
181  h->m[1] = p_Copy(q,r);
182#ifdef PDEBUG
183  Print("running syzygy comp. for nc_GCD:\n");
184#endif
185  ideal sh = idSyzygies(h, testHomog, &w);
186#ifdef PDEBUG
187  Print("done syzygy comp. for nc_GCD\n");
188#endif
189  /* in comm case, there is only 1 syzygy */
190  /*   singclap_gcd(); */
191  poly K, K1, K2;
192  K  = sh->m[0]; /* take just the first element - to be enhanced later */
193  K1 = pTakeOutComp(&K, 1); // 1st component is taken out from K
194//  pShift(&K,-2); // 2nd component to 0th comp.
195  K2 = pTakeOutComp(&K, 1);
196//  K2 = K;
197
198  Print("syz1: "); p_wrp(K1,r);
199  Print("syz2: "); p_wrp(K2,r);
200
201  /* checking signs before multiplying */   
202  number ck1 = p_GetCoeff(K1,r);
203  number ck2 = p_GetCoeff(K2,r);
204  BOOLEAN bck1, bck2;
205  bck1 = n_GreaterZero(ck1,r);
206  bck2 = n_GreaterZero(ck2,r);
207  /* K1 <0, K2 <0 (-K1,-K2)    */
208//   if ( !(bck1 && bck2) ) /* - , - */
209//   {
210//     K1 = p_Neg(K1,r);
211//     K2 = p_Neg(K2,r);
212//   }
213  id_Delete(&h,r);
214  h = idInit(2,1);
215  h->m[0] = p_Copy(K1,r);
216  h->m[1] = p_Copy(K2,r);
217  id_Delete(&sh,r);
218  return(h);
219}
220
221/* returns ideal (u,v) s.t. up + vq = 0 */
222
223ideal ncGCD(poly p, poly q, const ring r)
224{
225  // destroys p and q
226  // assume: p,q are in the comm. ring
227  // to be used in the coeff business
228#ifdef PDEBUG
229  PrintS(" GCD_start:");
230#endif
231  poly g = singclap_gcd(p_Copy(p,r),p_Copy(q,r));
232#ifdef PDEBUG
233  p_wrp(g,r);
234  PrintS(" GCD_end;\n");
235#endif
236  poly u = singclap_pdivide(q,g); //q/g
237  poly v = singclap_pdivide(p,g); //p/g
238  v = p_Neg(v,r);
239  p_Delete(&p,r);
240  p_Delete(&q,r);
241  ideal h = idInit(2,1);
242  h->m[0] = u; // p_Copy(u,r);
243  h->m[1] = v; // p_Copy(v,r);
244  return(h);
245}
246
247/* PINLINE1 void p_ExpVectorDiff
248   remains as is -> BUT we can do memory shift on smaller number of exp's */
249
250
251/*4 - follow the numbering of gring.cc
252* creates the S-polynomial of p1 and p2
253* do not destroy p1 and p2
254*/
255// poly nc_rat_CreateSpoly(poly p1, poly p2, poly spNoether, int ishift, const ring r)
256// {
257//   if ((p_GetComp(p1,r)!=p_GetComp(p2,r))
258//   && (p_GetComp(p1,r)!=0)
259//   && (p_GetComp(p2,r)!=0))
260//   {
261// #ifdef PDEBUG
262//     Print("nc_CreateSpoly : different components!");
263// #endif
264//     return(NULL);
265//   }
266//   /* prod. crit does not apply yet */
267// //   if ((r->nc->type==nc_lie) && pHasNotCF(p1,p2)) /* prod crit */
268// //   {
269// //     return(nc_p_Bracket_qq(pCopy(p2),p1));
270// //   }
271//   poly pL=pOne();
272//   poly m1=pOne();
273//   poly m2=pOne();
274//   /* define shift */
275//   int is = ishift; /* TODO */
276//   pLcmRat(p1,p2,pL,is);
277//   p_Setm(pL,r);
278//   poly pr1 = p_GetExp_k_n(p1,1,ishift-1,r); /* rat D-exp of p1 */
279//   poly pr2 = p_GetExp_k_n(p2,1,ishift-1,r); /* rat D-exp of p2 */
280// #ifdef PDEBUG
281//   p_Test(pL,r);
282// #endif
283//   p_ExpVectorDiff(m1,pL,p1,r); /* purely in D part by construction */
284//   //p_SetComp(m1,0,r);
285//   //p_Setm(m1,r);
286// #ifdef PDEBUG
287//   p_Test(m1,r);
288// #endif
289//   p_ExpVectorDiff(m2,pL,p2,r); /* purely in D part by construction */
290//   //p_SetComp(m2,0,r);
291//   //p_Setm(m2,r);
292// #ifdef PDEBUG
293//   p_Test(m2,r);
294// #endif
295//   p_Delete(&pL,r);
296//   /* zero exponents ! */
297
298//   /* EXTRACT LEADCOEF */
299
300//   poly H1  = p_HeadRat(p1,is,r);
301//   poly M1  = r->nc->p_Procs.mm_Mult_p(m1,p_Copy(H1,r),r);
302
303//   /* POLY:  number C1  = n_Copy(p_GetCoeff(M1,r),r); */
304//   /* RAT: */
305
306//   poly C1  = p_GetCoeffRat(M1,ishift,r);
307
308//   poly H2  = p_HeadRat(p2,is,r);
309//   poly M2  = r->nc->p_Procs.mm_Mult_p(m2,p_Copy(H2,r),r);
310
311//   /* POLY:  number C2  = n_Copy(p_GetCoeff(M2,r),r); */
312//   /* RAT: */
313
314//   poly C2  = p_GetCoeffRat(M2,ishift,r);
315
316// /* we do not assume that X's commute */
317// /* we just run NC syzygies */
318
319// /* NEW IDEA: change the ring to K<X>, map things there
320//    and return the result back; seems to be a good optimization */
321// /* to be done later */
322// /* problem: map to subalgebra. contexts, induced (non-unique) orderings etc. */
323
324//   intvec *w = NULL;
325//   ideal h = idInit(2,1);
326//   h->m[0] = p_Copy(C1,r);
327//   h->m[1] = p_Copy(C2,r);
328// #ifdef PDEBUG
329//   Print("running syzygy comp. for coeffs");
330// #endif
331//   ideal sh = idSyzygies(h, testHomog, &w);
332//   /* in comm case, there is only 1 syzygy */
333//   /*   singclap_gcd(); */
334//   poly K,K1,K2;
335//   K  = sh->m[0];
336//   K1 = pTakeOutComp(&K, 1); // 1st component is taken out from K
337//   pShift(&K,-2); // 2nd component to 0th comp.
338//   K2 = K;
339
340//   /* checking signs before multiplying */   
341//   number ck1 = p_GetCoeff(K1,r);
342//   number ck2 = p_GetCoeff(K2,r);
343//   BOOLEAN bck1, bck2;
344//   bck1 = n_GreaterZero(ck1,r);
345//   bck2 = n_GreaterZero(ck2,r);
346//   /* K1 >0, K2 >0 (K1,-K2)    */
347//   /* K1 >0, K2 <0 (K1,-K2)    */
348//   /* K1 <0, K2 >0 (-K1,K2)    */
349//   /* K1 <0, K2 <0 (-K1,K2)    */
350//   if ( (bck1) && (bck2) ) /* +, + */
351//   {
352//     K2 = p_Neg(K2,r);
353//   }
354//   if ( (bck1) && (!bck2) ) /* + , - */
355//   {
356//     K2 = p_Neg(K2,r);
357//   }
358//   if ( (!bck1) && (bck2) ) /* - , + */
359//   {
360//     K1 = p_Neg(K1,r);
361//   }
362//   if ( !(bck1 && bck2) ) /* - , - */
363//   {
364//     K1 = p_Neg(K1,r);
365//   }
366
367//   poly P1,P2;
368
369//   //  p_LmDeleteRat(M1,ishift,r); // get tail(D^(gamma-alpha) * lm(p1)) = h_f
370//   P1 = p_Copy(p1,r);
371//   p_LmDeleteAndNextRat(P1,ishift,r); // get tail(p1) = t_f
372//   P1 = r->nc->p_Procs.mm_Mult_p(m1,P1,r);
373//   P1 = p_Add_q(P1,M1,r);
374
375//   //  p_LmDeleteRat(M2,ishift,r);
376//   P2 = p_Copy(p2,r);
377//   p_LmDeleteAndNextRat(P2,ishift,r);// get tail(p2)=t_g
378//   P2 = r->nc->p_Procs.mm_Mult_p(m2,P2,r);
379//   P2 = p_Add_q(P2,M2,r);
380
381//   /* coeff business */
382
383//   P1 = p_Mult_q(P1,K1,r);
384//   P2 = p_Mult_q(P2,K2,r);
385//   P1 = p_Add_q(P1,P2,r);
386
387//   /* cleaning up */
388
389// #ifdef PDEBUG
390//   p_Test(p1,r);
391// #endif
392//   /* questionable: */
393//   if (P1!=NULL) pCleardenom(P1);
394//   if (P1!=NULL) pContent(P1);
395//   return(P1);
396// }
397
398
399/*4 - follow the numbering of gring.cc
400* creates the S-polynomial of p1 and p2
401* do not destroy p1 and p2
402*/
403poly nc_rat_CreateSpoly(poly pp1, poly pp2, int ishift, const ring r)
404{
405
406  poly p1 = p_Copy(pp1,r);
407  poly p2 = p_Copy(pp2,r);
408
409  const long lCompP1 = p_GetComp(p1,r);
410  const long lCompP2 = p_GetComp(p2,r);
411
412  if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
413  {
414#ifdef PDEBUG
415    Werror("nc_rat_CreateSpoly: different non-zero components!");
416#endif
417    return(NULL);
418  }
419
420/* note: prod. crit does not apply! */
421  poly pL=pOne();
422  poly m1=pOne();
423  poly m2=pOne();
424  int is = ishift; /* TODO */
425  pLcmRat(p1,p2,pL,is);
426  p_Setm(pL,r);
427#ifdef PDEBUG
428  p_Test(pL,r);
429#endif
430  poly pr1 = p_GetExp_k_n(p1,1,ishift,r); /* rat D-exp of p1 */
431  poly pr2 = p_GetExp_k_n(p2,1,ishift,r); /* rat D-exp of p2 */
432  p_ExpVectorDiff(m1,pL,pr1,r); /* purely in D part by construction */
433  p_ExpVectorDiff(m2,pL,pr2,r); /* purely in D part by construction */
434  p_Delete(&pr1,r);
435  p_Delete(&pr2,r);
436  p_Delete(&pL,r);
437#ifdef PDEBUG
438  p_Test(m1,r);
439  PrintS("d^{gamma-alpha} = "); p_wrp(m1,r); PrintLn();
440  p_Test(m2,r);
441  PrintS("d^{gamma-beta} = "); p_wrp(m2,r); PrintLn();
442#endif
443
444  poly HF = NULL;
445  HF = p_HeadRat(p1,is,r); // lm_D(f)
446  HF  = nc_mm_Mult_p(m1, HF, r); // // d^{gamma-alpha} lm_D(f)
447  poly C  = p_GetCoeffRat(HF,  is, r); // c = lc_D(h_f) in the paper
448
449  poly HG = NULL;
450  HG = p_HeadRat(p2,is,r); // lm_D(g)
451  HG  = nc_mm_Mult_p(m2, HG, r); // // d^{gamma-beta} lm_D(g)
452  poly K  = p_GetCoeffRat(HG,  is, r); // k = lc_D(h_g) in the paper
453
454#ifdef PDEBUG
455  PrintS("f: "); p_wrp(p1,r); PrintS("\n");
456  PrintS("c: "); p_wrp(C,r); PrintS("\n");
457  PrintS("g: "); p_wrp(p2,r); PrintS("\n");
458  PrintS("k: "); p_wrp(K,r); PrintS("\n");
459#endif
460 
461  ideal ncsyz = ncGCD(C,K,r);
462  poly KK = ncsyz->m[0]; ncsyz->m[0]=NULL; //p_Copy(ncsyz->m[0],r); // k'
463  poly CC = ncsyz->m[1]; ncsyz->m[1]= NULL; //p_Copy(ncsyz->m[1],r); // c'
464  id_Delete(&ncsyz,r);
465
466  p_LmDeleteAndNextRat(&p1, is+1, r); // t_f
467  p_LmDeleteAndNextRat(&HF, is+1, r); // r_f = h_f - lt_D(h_f)
468
469  p_LmDeleteAndNextRat(&p2, is+1, r); // t_g
470  p_LmDeleteAndNextRat(&HG, is+1, r); // r_g = h_g - lt_D(h_g)
471
472
473#ifdef PDEBUG
474  PrintS(" t_f: "); p_wrp(p1,r); PrintS("\n"); 
475  PrintS(" t_g: "); p_wrp(p2,r); PrintS("\n"); 
476  PrintS(" r_f: "); p_wrp(HF,r); PrintS("\n"); 
477  PrintS(" r_g: "); p_wrp(HG,r); PrintS("\n"); 
478  PrintS(" c': "); p_wrp(CC,r); PrintS("\n"); 
479  PrintS(" k': "); p_wrp(KK,r); PrintS("\n"); 
480
481#endif
482
483  // k'(r_f + d^{gamma-alpha} t_f)
484
485  p1 = p_Mult_q(m1, p1, r); // p1 = d^{gamma-alpha} t_f
486  p1 = p_Add_q(p1,HF,r); // p1 = r_f + d^{gamma-alpha} t_f
487  p1 = p_Mult_q(KK,p1,r); // p1 = k'(r_f + d^{gamma-alpha} t_f)
488
489  // c'(r_f + d^{gamma-beta} t_g)
490
491  p2 = p_Mult_q(m2, p2, r); // p2 = d^{gamma-beta} t_g
492  p2 = p_Add_q(p2,HG,r); // p2 = r_g + d^{gamma-beta} t_g
493  p2 = p_Mult_q(CC,p2,r); // p2 = c'(r_g + d^{gamma-beta} t_g)
494
495#ifdef PDEBUG
496  p_Test(p1,r);
497  p_Test(p2,r);
498  PrintS(" k'(r_f + d^{gamma-alpha} t_f): "); p_wrp(p1,r);
499  PrintS(" c'(r_g + d^{gamma-beta} t_g): "); p_wrp(p2,r);
500#endif
501
502  poly out = p_Add_q(p1,p2,r); // delete p1, p2; // the sum
503
504#ifdef PDEBUG
505  p_Test(out,r);
506#endif
507
508  if ( out!=NULL ) pContent(out);
509  return(out);
510}
511
512
513/*2
514* reduction of p2 with p1
515* do not destroy p1, but p2
516* p1 divides p2 -> for use in NF algorithm
517* works in an integer fashion
518*/
519
520poly nc_rat_ReduceSpolyNew(const poly p1, poly p2, int ishift, const ring r)
521{
522  const long lCompP1 = p_GetComp(p1,r);
523  const long lCompP2 = p_GetComp(p2,r);
524
525  if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
526  {
527#ifdef PDEBUG
528    Werror("nc_rat_ReduceSpolyNew: different non-zero components!");
529#endif
530    return(NULL);
531  }
532
533  int is = ishift; /* TODO */
534
535  poly m = pOne();
536  p_ExpVectorDiffRat(m, p2, p1, ishift, r); // includes X and D parts
537  //p_Setm(m,r);
538  //  m = p_GetExp_k_n(m,1,ishift,r); /* rat D-exp of m */
539#ifdef PDEBUG
540  p_Test(m,r);
541  PrintS("d^alpha = "); p_wrp(m,r); PrintLn();
542#endif
543
544  /* pSetComp(m,r)=0? */
545  poly HH = NULL;
546  poly H  = NULL;
547  HH = p_HeadRat(p1,is,r); //p_Copy(p_HeadRat(p1,is,r),r); // lm_D(g)
548//  H  = r->nc->p_Procs.mm_Mult_p(m, p_Copy(HH, r), r); // d^aplha lm_D(g)
549  H  = nc_mm_Mult_p(m, HH, r); // d^aplha lm_D(g) == h_g in the paper
550
551  poly K  = p_GetCoeffRat(H,  is, r); //p_Copy( p_GetCoeffRat(H,  is, r), r); // k in the paper
552  poly P  = p_GetCoeffRat(p2, is, r); //p_Copy( p_GetCoeffRat(p2, is, r), r); // lc_D(p_2) == lc_D(f)
553
554#ifdef PDEBUG
555  PrintS("k: "); p_wrp(K,r); PrintS("\n");
556  PrintS("p: "); p_wrp(P,r); PrintS("\n");
557  PrintS("f: "); p_wrp(p2,r); PrintS("\n");
558  PrintS("g: "); p_wrp(p1,r); PrintS("\n");
559#endif
560  // alt:
561  poly out = p_Copy(p1,r);
562  p_LmDeleteAndNextRat(&out, is+1, r); // out == t_g
563
564  ideal ncsyz = ncGCD(P,K,r);
565  poly KK = ncsyz->m[0]; ncsyz->m[0]=NULL; //p_Copy(ncsyz->m[0],r); // k'
566  poly PP = ncsyz->m[1]; ncsyz->m[1]= NULL; //p_Copy(ncsyz->m[1],r); // p'
567
568#ifdef PDEBUG
569  PrintS("t_g: "); p_wrp(out,r);
570  PrintS("k': "); p_wrp(KK,r); PrintS("\n"); 
571  PrintS("p': "); p_wrp(PP,r); PrintS("\n"); 
572#endif
573  id_Delete(&ncsyz,r);
574  p_LmDeleteAndNextRat(&p2, is+1, r); // t_f
575  p_LmDeleteAndNextRat(&H, is+1, r); // r_g = h_g - lt_D(h_g)
576
577#ifdef PDEBUG
578  PrintS(" t_f: "); p_wrp(p2,r);
579  PrintS(" r_g: "); p_wrp(H,r);
580#endif
581
582  p2 = p_Mult_q(KK, p2, r); // p2 = k' t_f
583
584#ifdef PDEBUG
585  p_Test(p2,r);
586  PrintS(" k' t_f: "); p_wrp(p2,r);
587#endif
588
589//  out = r->nc->p_Procs.mm_Mult_p(m, out, r); // d^aplha t_g
590  out = nc_mm_Mult_p(m, out, r); // d^aplha t_g 
591  p_Delete(&m,r);
592
593#ifdef PDEBUG
594  PrintS(" d^a t_g: "); p_wrp(out,r);
595  PrintS(" end reduction\n");
596#endif
597
598  out = p_Add_q(H, out, r); // r_g + d^a t_g
599
600#ifdef PDEBUG
601  p_Test(out,r);
602#endif
603  out = p_Mult_q(PP, out, r); // p' (r_g + d^a t_g)
604  out = p_Add_q(p2,out,r); // delete out, p2; // the sum
605
606#ifdef PDEBUG
607  p_Test(out,r);
608#endif
609
610  if ( out!=NULL ) pContent(out);
611  return(out);
612}
613
614// return: FALSE, if there exists i in ishift..r->N,
615//                 such that a->exp[i] > b->exp[i]
616//         TRUE, otherwise
617
618BOOLEAN p_DivisibleByRat(poly a, poly b, int ishift, const ring r)
619{
620#ifdef PDEBUG
621  PrintS("invoke p_DivByRat with a = ");
622  p_wrp(p_Head(a,r),r);
623  PrintS(" and b= ");
624  p_wrp(p_Head(b,r),r); 
625  PrintLn();
626#endif
627  int i;
628  for(i=r->N; i>ishift; i--)
629  {
630#ifdef PDEBUG
631    Print("i=%d,",i);
632#endif
633    if (p_GetExp(a,i,r) > p_GetExp(b,i,r)) return FALSE;
634  }
635  return ((p_GetComp(a,r)==p_GetComp(b,r)) || (p_GetComp(a,r)==0));
636}
637/*2
638*reduces h with elements from reducer choosing the best possible
639* element in t with respect to the given red_length
640* arrays reducer and red_length are [0..(rl-1)]
641*/
642int redRat (poly* h, poly *reducer, int *red_length, int rl, int ishift, ring r)
643{
644  if ((*h)==NULL) return 0;
645
646  int j,i,l;
647
648  loop
649  {
650    j=rl;l=MAX_INT_VAL;
651    for(i=rl-1;i>=0;i--)
652    {
653      //      Print("test %d, l=%d (curr=%d, l=%d\n",i,red_length[i],j,l);
654      if ((l>red_length[i]) && (p_DivisibleByRat(reducer[i],*h,ishift,r)))
655      {
656        j=i; l=red_length[i];
657        //        PrintS(" yes\n");
658      }
659      //      else PrintS(" no\n");
660    }
661    if (j >=rl)
662    {
663      return 1; // not reducible
664    }
665
666    if (TEST_OPT_DEBUG)
667    {
668      PrintS("reduce ");
669      p_wrp(*h,r);
670      PrintS(" with ");
671      p_wrp(reducer[j],r);
672    }
673    poly hh=nc_rat_ReduceSpolyNew(reducer[j], *h, ishift, r);
674    //    p_Delete(h,r);
675    *h=hh;
676    if (TEST_OPT_DEBUG)
677    {
678      PrintS(" to ");
679      p_wrp(*h,r);
680      PrintLn();
681    }
682    if ((*h)==NULL)
683    {
684      return 0;
685    }
686  }
687}
688
689poly pContentRat(poly ph)
690// destroys ph
691// for rat coefficients in K(x1,..xN)
692{
693
694  // init array of RatLeadCoeffs
695  //  poly p_GetCoeffRat(poly p, int ishift, ring r);
696
697  int len=pLength(ph);
698  poly *C = (poly *)omAlloc0((len+1)*sizeof(poly));  //rat coeffs
699  poly *LM = (poly *)omAlloc0((len+1)*sizeof(poly));  // rat lead terms
700  int *D = (int *)omAlloc0((len+1)*sizeof(int));  //degrees of coeffs
701  int *L = (int *)omAlloc0((len+1)*sizeof(int));  //lengths of coeffs
702  int k = 0;
703  poly p = ph;
704  int HasConstantCoef = 0;
705  int mintdeg = pTotaldegree(p);
706  int minlen = len;
707  while (p!=NULL)
708  {
709    LM[k] = p_HeadRat(p, currRing->real_var_start, currRing); 
710    C[k] = p_GetCoeffRat(p, currRing->real_var_start, currRing);
711    D[k] =  pTotaldegree(C[k]);
712    mintdeg = si_min(mintdeg,D[k]);
713    L[k] = pLength(C[k]);
714    minlen = si_min(minlen,L[k]);
715    if (pIsConstant(C[k]))
716    {
717      // C[k] = const, so the content will be numerical
718      HasConstantCoef = 1;
719      // smth like goto cleanup and return(pContent(p));
720    }
721    p_LmDeleteAndNextRat(&p, currRing->real_var_start, currRing);
722    k++;
723  }
724
725  // look for 1 element of minimal degree and of minimal length
726  k--;
727  int dd = 0; int i;
728  int mindeglen = len;
729  int pmindeglen;
730  for(i=0; i<=k; i++)
731  {
732    if (D[i] == mintdeg)
733    {
734      if (L[i] < mindeglen)
735      {
736        mindeglen=L[i];
737        pmindeglen = i;
738      }
739    }
740  }
741  poly d = pCopy(C[pmindeglen]);
742  // there are dd>=1 mindeg elements
743  // and pmideglen is the coordinate of one of the smallest among them
744
745  //  poly g = singclap_gcd(p_Copy(p,r),p_Copy(q,r));
746  //  return naGcd(d,d2,currRing);
747
748  // adjoin pContentRat here?
749  for(i=0; i<=k; i++)
750  {
751    d=singclap_gcd(d,pCopy(C[i]));
752    if (pTotaldegree(d)==0) 
753    {
754      // cleanup, pContent, return
755      pDelete(&d);
756      for(;k>=0;k--)
757      {
758        pDelete(&C[k]);
759        pDelete(&LM[k]);
760      }
761      p=pContent(ph);
762      goto cleanup;
763    }
764  }
765  for(i=0; i<=k; i++)
766  {
767   poly h=singclap_pdivide(C[i],d);
768   pDelete(&C[i]);
769   C[i]=h;
770  }
771 
772  // zusammensetzen,
773  p=NULL; // just to be sure
774  for(i=0; i<=k; i++)
775  {
776   p=pAdd(p,pMult(C[i],LM[i]));
777   C[i]=NULL; LM[i]=NULL;
778  }
779  // aufraeumen, return
780cleanup:
781  omFree(C);
782  omFree(LM);
783  omFree(D);
784  omFree(L);
785  return p;
786}
787#endif
Note: See TracBrowser for help on using the repository browser.