source: git/kernel/ratgring.cc @ 341696

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