source: git/kernel/ratgring.cc @ ad36a1

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