source: git/kernel/ratgring.cc @ a82c308

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