source: git/kernel/ratgring.cc @ fbc7cb

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