source: git/kernel/ratgring.cc @ 019a78

spielwiese
Last change on this file since 019a78 was 019a78, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: syntax git-svn-id: file:///usr/local/Singular/svn/trunk@11427 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    ratgring.cc
6 *  Purpose: Ore-noncommutative kernel procedures
7 *  Author:  levandov (Viktor Levandovsky)
8 *  Created: 8/00 - 11/00
9 *  Version: $Id: ratgring.cc,v 1.13 2009-02-21 16:01:05 Singular Exp $
10 *******************************************************************/
11#include "mod2.h"
12#include "ratgring.h"
13#ifdef HAVE_RATGRING
14#include "gring.h"
15#include "febase.h"
16#include "ring.h"
17#include "polys.h"
18#include "numbers.h"
19#include "ideals.h"
20#include "matpol.h"
21#include "kbuckets.h"
22#include "kstd1.h"
23#include "sbuckets.h"
24#include "prCopy.h"
25#include "p_Mult_q.h"
26#include "clapsing.h"
27
28void pLcmRat(poly a, poly b, poly m, int rat_shift)
29{
30  /* rat_shift is the last exp one should count with */
31  int i;
32  for (i=pVariables; i>=rat_shift; i--)
33  {
34    pSetExp(m,i, si_max( pGetExp(a,i), pGetExp(b,i)));
35  }
36  pSetComp(m, si_max(pGetComp(a), pGetComp(b)));
37  /* Don't do a pSetm here, otherwise hres/lres chockes */ 
38}
39
40// void pLcmRat(poly a, poly b, poly m, poly pshift)
41// {
42//   /* shift is the exp of rational elements */
43//   int i;
44//   for (i=pVariables; i; i--)
45//   {
46//     if (!pGetExp(pshift,i))
47//     {
48//       pSetExp(m,i, si_max( pGetExp(a,i), pGetExp(b,i)));
49//     }
50//     else
51//     {
52//       /* do we really need it? */
53//       pSetExp(m,i,0);
54//     }
55//   }
56//   pSetComp(m, si_max(pGetComp(a), pGetComp(b)));
57//   /* Don't do a pSetm here, otherwise hres/lres chockes */
58// }
59
60/* returns a subpoly of p, s.t. its monomials have the same D-part */
61
62poly p_HeadRat(poly p, int ishift, ring r)
63{
64  poly q   = pNext(p);
65  if (q == NULL) return p;
66  poly res = p_Head(p,r);
67  while ( (q!=NULL) && (p_Comp_k_n(p, q, ishift+1, r)))
68  {
69    res = p_Add_q(res,p_Head(q,r),r);
70    q   = pNext(q);
71  }
72  return res;
73}
74
75/* TO TEST!!! */
76/* returns x-coeff of p, i.e. a poly in x, s.t. corresponding xd-monomials
77have the same D-part */
78
79poly p_GetCoeffRat(poly p, int ishift, ring r)
80{
81  poly q   = pNext(p);
82  poly res; // = p_Head(p,r);
83  res = p_GetExp_k_n(p, ishift+1, r->N, r);
84  p_SetCoeff(res,n_Copy(p_GetCoeff(p,r),r),r);
85  poly s;
86  while ((q!= NULL) && (p_Comp_k_n(p, q, ishift+1, r)))
87  {
88    s   = p_GetExp_k_n(q, ishift+1, r->N, r);
89    p_SetCoeff(s,n_Copy(p_GetCoeff(q,r),r),r);
90    res = p_Add_q(res,s,r);
91    q   = pNext(q);
92  }
93  return res;
94}
95
96void p_LmDeleteAndNextRat(poly *p, int ishift, ring r)
97{
98  /* modifies p*/
99  //  Print("start: "); Print(" "); p_wrp(*p,r);
100  p_LmCheckPolyRing2(*p, r);
101  poly q = p_Head(*p,r);
102  // in the next line ishift is correct
103  while ( ( (*p)!=NULL ) && ( p_Comp_k_n(*p, q, ishift, r) ))
104  {
105    p_LmDelete(p,r);
106    //    Print("while: ");p_wrp(*p,r);Print(" ");
107  }
108  //  p_wrp(*p,r);Print(" ");
109  //  PrintS("end\n");
110  p_LmDelete(&q,r);
111}
112
113/* to test!!! */
114/* ExpVector(pr) = ExpVector(p1) - ExpVector(p2) */
115void p_ExpVectorDiffRat(poly pr, poly p1, poly p2, int ishift, ring r)
116{
117  p_LmCheckPolyRing1(p1, r);
118  p_LmCheckPolyRing1(p2, r);
119  p_LmCheckPolyRing1(pr, r); 
120  int i;
121  poly t=pr;
122  Exponent_t e1,e2;
123  for (i=ishift+1; i<=r->N; i++)
124  {
125    e1 = p_GetExp(p1, i, r);
126    e2 = p_GetExp(p2, i, r);
127    //    pAssume1(p_GetExp(p1, i, r) >= p_GetExp(p2, i, r));
128    if (e1 < e2)
129    {
130#ifdef PDEBUG
131      Print("negative ExpVectorDiff\n");
132#endif   
133      p_Delete(&t,r);
134      break;
135    }
136    else
137    {
138      p_SetExp(t,i, e1-e2,r);
139    }
140  }
141  p_Setm(t,r);
142}
143
144/* returns ideal (u,v) s.t. up + vq = 0 */
145
146ideal ncGCD2(poly p, poly q, const ring r)
147{
148  // todo: must destroy p,q
149  intvec *w = NULL;
150  ideal h = idInit(2,1);
151  h->m[0] = p_Copy(p,r);
152  h->m[1] = p_Copy(q,r);
153#ifdef PDEBUG
154  Print("running syzygy comp. for nc_GCD:\n");
155#endif
156  ideal sh = idSyzygies(h, testHomog, &w);
157#ifdef PDEBUG
158  Print("done syzygy comp. for nc_GCD\n");
159#endif
160  /* in comm case, there is only 1 syzygy */
161  /*   singclap_gcd(); */
162  poly K, K1, K2;
163  K  = sh->m[0]; /* take just the first element - to be enhanced later */
164  K1 = pTakeOutComp(&K, 1); // 1st component is taken out from K
165//  pShift(&K,-2); // 2nd component to 0th comp.
166  K2 = pTakeOutComp(&K, 1);
167//  K2 = K;
168
169  Print("syz1: "); p_wrp(K1,r);
170  Print("syz2: "); p_wrp(K2,r);
171
172  /* checking signs before multiplying */   
173  number ck1 = p_GetCoeff(K1,r);
174  number ck2 = p_GetCoeff(K2,r);
175  BOOLEAN bck1, bck2;
176  bck1 = n_GreaterZero(ck1,r);
177  bck2 = n_GreaterZero(ck2,r);
178  /* K1 <0, K2 <0 (-K1,-K2)    */
179//   if ( !(bck1 && bck2) ) /* - , - */
180//   {
181//     K1 = p_Neg(K1,r);
182//     K2 = p_Neg(K2,r);
183//   }
184  id_Delete(&h,r);
185  h = idInit(2,1);
186  h->m[0] = p_Copy(K1,r);
187  h->m[1] = p_Copy(K2,r);
188  id_Delete(&sh,r);
189  return(h);
190}
191
192/* returns ideal (u,v) s.t. up + vq = 0 */
193
194ideal ncGCD(poly p, poly q, const ring r)
195{
196  // destroys p and q
197  // assume: p,q are in the comm. ring
198  // to be used in the coeff business
199#ifdef PDEBUG
200  PrintS(" GCD_start:");
201#endif
202  poly g = singclap_gcd(p_Copy(p,r),p_Copy(q,r));
203#ifdef PDEBUG
204  p_wrp(g,r);
205  PrintS(" GCD_end;\n");
206#endif
207  poly u = singclap_pdivide(q,g); //q/g
208  poly v = singclap_pdivide(p,g); //p/g
209  p_Delete(&p,r);
210  p_Delete(&q,r);
211  ideal h = idInit(2,1);
212  h->m[0] = u; // p_Copy(u,r);
213  h->m[1] = v; // p_Copy(v,r);
214  return(h);
215}
216
217/* PINLINE1 void p_ExpVectorDiff
218   remains as is -> BUT we can do memory shift on smaller number of exp's */
219
220
221/*4 - follow the numbering of gring.cc
222* creates the S-polynomial of p1 and p2
223* do not destroy p1 and p2
224*/
225// poly nc_rat_CreateSpoly(poly p1, poly p2, poly spNoether, int ishift, const ring r)
226// {
227//   if ((p_GetComp(p1,r)!=p_GetComp(p2,r))
228//   && (p_GetComp(p1,r)!=0)
229//   && (p_GetComp(p2,r)!=0))
230//   {
231// #ifdef PDEBUG
232//     Print("nc_CreateSpoly : different components!");
233// #endif
234//     return(NULL);
235//   }
236//   /* prod. crit does not apply yet */
237// //   if ((r->nc->type==nc_lie) && pHasNotCF(p1,p2)) /* prod crit */
238// //   {
239// //     return(nc_p_Bracket_qq(pCopy(p2),p1));
240// //   }
241//   poly pL=pOne();
242//   poly m1=pOne();
243//   poly m2=pOne();
244//   /* define shift */
245//   int is = ishift; /* TODO */
246//   pLcmRat(p1,p2,pL,is);
247//   p_Setm(pL,r);
248//   poly pr1 = p_GetExp_k_n(p1,1,ishift-1,r); /* rat D-exp of p1 */
249//   poly pr2 = p_GetExp_k_n(p2,1,ishift-1,r); /* rat D-exp of p2 */
250// #ifdef PDEBUG
251//   p_Test(pL,r);
252// #endif
253//   p_ExpVectorDiff(m1,pL,p1,r); /* purely in D part by construction */
254//   //p_SetComp(m1,0,r);
255//   //p_Setm(m1,r);
256// #ifdef PDEBUG
257//   p_Test(m1,r);
258// #endif
259//   p_ExpVectorDiff(m2,pL,p2,r); /* purely in D part by construction */
260//   //p_SetComp(m2,0,r);
261//   //p_Setm(m2,r);
262// #ifdef PDEBUG
263//   p_Test(m2,r);
264// #endif
265//   p_Delete(&pL,r);
266//   /* zero exponents ! */
267
268//   /* EXTRACT LEADCOEF */
269
270//   poly H1  = p_HeadRat(p1,is,r);
271//   poly M1  = r->nc->p_Procs.mm_Mult_p(m1,p_Copy(H1,r),r);
272
273//   /* POLY:  number C1  = n_Copy(p_GetCoeff(M1,r),r); */
274//   /* RAT: */
275
276//   poly C1  = p_GetCoeffRat(M1,ishift,r);
277
278//   poly H2  = p_HeadRat(p2,is,r);
279//   poly M2  = r->nc->p_Procs.mm_Mult_p(m2,p_Copy(H2,r),r);
280
281//   /* POLY:  number C2  = n_Copy(p_GetCoeff(M2,r),r); */
282//   /* RAT: */
283
284//   poly C2  = p_GetCoeffRat(M2,ishift,r);
285
286// /* we do not assume that X's commute */
287// /* we just run NC syzygies */
288
289// /* NEW IDEA: change the ring to K<X>, map things there
290//    and return the result back; seems to be a good optimization */
291// /* to be done later */
292// /* problem: map to subalgebra. contexts, induced (non-unique) orderings etc. */
293
294//   intvec *w = NULL;
295//   ideal h = idInit(2,1);
296//   h->m[0] = p_Copy(C1,r);
297//   h->m[1] = p_Copy(C2,r);
298// #ifdef PDEBUG
299//   Print("running syzygy comp. for coeffs");
300// #endif
301//   ideal sh = idSyzygies(h, testHomog, &w);
302//   /* in comm case, there is only 1 syzygy */
303//   /*   singclap_gcd(); */
304//   poly K,K1,K2;
305//   K  = sh->m[0];
306//   K1 = pTakeOutComp(&K, 1); // 1st component is taken out from K
307//   pShift(&K,-2); // 2nd component to 0th comp.
308//   K2 = K;
309
310//   /* checking signs before multiplying */   
311//   number ck1 = p_GetCoeff(K1,r);
312//   number ck2 = p_GetCoeff(K2,r);
313//   BOOLEAN bck1, bck2;
314//   bck1 = n_GreaterZero(ck1,r);
315//   bck2 = n_GreaterZero(ck2,r);
316//   /* K1 >0, K2 >0 (K1,-K2)    */
317//   /* K1 >0, K2 <0 (K1,-K2)    */
318//   /* K1 <0, K2 >0 (-K1,K2)    */
319//   /* K1 <0, K2 <0 (-K1,K2)    */
320//   if ( (bck1) && (bck2) ) /* +, + */
321//   {
322//     K2 = p_Neg(K2,r);
323//   }
324//   if ( (bck1) && (!bck2) ) /* + , - */
325//   {
326//     K2 = p_Neg(K2,r);
327//   }
328//   if ( (!bck1) && (bck2) ) /* - , + */
329//   {
330//     K1 = p_Neg(K1,r);
331//   }
332//   if ( !(bck1 && bck2) ) /* - , - */
333//   {
334//     K1 = p_Neg(K1,r);
335//   }
336
337//   poly P1,P2;
338
339//   //  p_LmDeleteRat(M1,ishift,r); // get tail(D^(gamma-alpha) * lm(p1)) = h_f
340//   P1 = p_Copy(p1,r);
341//   p_LmDeleteAndNextRat(P1,ishift,r); // get tail(p1) = t_f
342//   P1 = r->nc->p_Procs.mm_Mult_p(m1,P1,r);
343//   P1 = p_Add_q(P1,M1,r);
344
345//   //  p_LmDeleteRat(M2,ishift,r);
346//   P2 = p_Copy(p2,r);
347//   p_LmDeleteAndNextRat(P2,ishift,r);// get tail(p2)=t_g
348//   P2 = r->nc->p_Procs.mm_Mult_p(m2,P2,r);
349//   P2 = p_Add_q(P2,M2,r);
350
351//   /* coeff business */
352
353//   P1 = p_Mult_q(P1,K1,r);
354//   P2 = p_Mult_q(P2,K2,r);
355//   P1 = p_Add_q(P1,P2,r);
356
357//   /* cleaning up */
358
359// #ifdef PDEBUG
360//   p_Test(p1,r);
361// #endif
362//   /* questionable: */
363//   if (P1!=NULL) pCleardenom(P1);
364//   if (P1!=NULL) pContent(P1);
365//   return(P1);
366// }
367
368/*2
369* reduction of p2 with p1
370* do not destroy p1, but p2
371* p1 divides p2 -> for use in NF algorithm
372* works in an integer fashion
373*/
374
375poly nc_rat_ReduceSpolyNew(const poly p1, poly p2, int ishift, const ring r)
376{
377  const long lCompP1 = p_GetComp(p1,r);
378  const long lCompP2 = p_GetComp(p2,r);
379
380  if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
381  {
382#ifdef PDEBUG
383    Werror("nc_rat_ReduceSpolyNew: different non-zero components!");
384#endif
385    return(NULL);
386  }
387
388  int is = ishift; /* TODO */
389
390  poly m = pOne();
391  p_ExpVectorDiffRat(m, p2, p1, ishift, r); // includes X and D parts
392  //p_Setm(m,r);
393  //  m = p_GetExp_k_n(m,1,ishift,r); /* rat D-exp of m */
394#ifdef PDEBUG
395  p_Test(m,r);
396  PrintS("d^alpha = "); p_wrp(m,r); PrintLn();
397#endif
398
399  /* pSetComp(m,r)=0? */
400  poly HH = NULL;
401  poly H  = NULL;
402  HH = p_HeadRat(p1,is,r); //p_Copy(p_HeadRat(p1,is,r),r); // lm_D(g)
403//  H  = r->nc->p_Procs.mm_Mult_p(m, p_Copy(HH, r), r); // d^aplha lm_D(g)
404  H  = nc_mm_Mult_p(m, HH, r); // d^aplha lm_D(g) == h_g in the paper
405
406  poly K  = p_GetCoeffRat(H,  is, r); //p_Copy( p_GetCoeffRat(H,  is, r), r); // k in the paper
407  poly P  = p_GetCoeffRat(p2, is, r); //p_Copy( p_GetCoeffRat(p2, is, r), r); // lc_D(p_2) == lc_D(f)
408
409#ifdef PDEBUG
410  PrintS("k: "); p_wrp(K,r); PrintS("\n");
411  PrintS("p: "); p_wrp(P,r); PrintS("\n");
412  PrintS("f: "); p_wrp(p2,r); PrintS("\n");
413  PrintS("g: "); p_wrp(p1,r); PrintS("\n");
414#endif
415  // alt:
416  poly out = p_Copy(p1,r);
417  p_LmDeleteAndNextRat(&out, is+1, r); // out == t_g
418
419  ideal ncsyz = ncGCD(P,K,r);
420  poly KK = ncsyz->m[0]; ncsyz->m[0]=NULL; //p_Copy(ncsyz->m[0],r); // k'
421  poly PP = ncsyz->m[1]; ncsyz->m[1]= NULL; //p_Copy(ncsyz->m[1],r); // p'
422
423#ifdef PDEBUG
424  PrintS("t_g: "); p_wrp(out,r);
425  PrintS("k': "); p_wrp(KK,r); PrintS("\n"); 
426  PrintS("p': "); p_wrp(PP,r); PrintS("\n"); 
427#endif
428  id_Delete(&ncsyz,r);
429  p_LmDeleteAndNextRat(&p2, is+1, r); // t_f
430  p_LmDeleteAndNextRat(&H, is+1, r); // r_g = h_g - lt_D(h_g)
431
432#ifdef PDEBUG
433  PrintS(" t_f: "); p_wrp(p2,r);
434  PrintS(" r_g: "); p_wrp(H,r);
435#endif
436
437  p2 = p_Mult_q(KK, p2, r); // p2 = k' t_f
438
439#ifdef PDEBUG
440  p_Test(p2,r);
441  PrintS(" k' t_f: "); p_wrp(p2,r);
442#endif
443
444//  out = r->nc->p_Procs.mm_Mult_p(m, out, r); // d^aplha t_g
445  out = nc_mm_Mult_p(m, out, r); // d^aplha t_g 
446  p_Delete(&m,r);
447
448#ifdef PDEBUG
449  PrintS(" d^a t_g: "); p_wrp(out,r);
450  PrintS(" end reduction\n");
451#endif
452
453  out = p_Add_q(H, out, r); // r_g + d^a t_g
454
455#ifdef PDEBUG
456  p_Test(out,r);
457#endif
458  out = p_Mult_q(PP, out, r); // p' (r_g + d^a t_g)
459  out = p_Add_q(p2,out,r); // delete out, p2; // the sum
460
461#ifdef PDEBUG
462  p_Test(out,r);
463#endif
464
465  if ( out!=NULL ) pContent(out);
466  return(out);
467}
468
469// return: FALSE, if there exists i in ishift..r->N,
470//                 such that a->exp[i] > b->exp[i]
471//         TRUE, otherwise
472
473BOOLEAN p_DivisibleByRat(poly a, poly b, int ishift, const ring r)
474{
475#ifdef PDEBUG
476  PrintS("invoke p_DivByRat with a = ");
477  p_wrp(p_Head(a,r),r);
478  PrintS(" and b= ");
479  p_wrp(p_Head(b,r),r); 
480  PrintLn();
481#endif
482  int i;
483  for(i=r->N; i>ishift; i--)
484  {
485#ifdef PDEBUG
486    Print("i=%d,",i);
487#endif
488    if (p_GetExp(a,i,r) > p_GetExp(b,i,r)) return FALSE;
489  }
490  return ((p_GetComp(a,r)==p_GetComp(b,r)) || (p_GetComp(a,r)==0));
491}
492/*2
493*reduces h with elements from reducer choosing the best possible
494* element in t with respect to the given red_length
495* arrays reducer and red_length are [0..(rl-1)]
496*/
497int redRat (poly* h, poly *reducer, int *red_length, int rl, int ishift, ring r)
498{
499  if ((*h)==NULL) return 0;
500
501  int j,i,l;
502
503  loop
504  {
505    j=rl;l=MAX_INT_VAL;
506    for(i=rl-1;i>=0;i--)
507    {
508      //      Print("test %d, l=%d (curr=%d, l=%d\n",i,red_length[i],j,l);
509      if ((l>red_length[i]) && (p_DivisibleByRat(reducer[i],*h,ishift,r)))
510      {
511        j=i; l=red_length[i];
512        //        PrintS(" yes\n");
513      }
514      //      else PrintS(" no\n");
515    }
516    if (j >=rl)
517    {
518      return 1; // not reducible
519    }
520
521    if (TEST_OPT_DEBUG)
522    {
523      PrintS("reduce ");
524      p_wrp(*h,r);
525      PrintS(" with ");
526      p_wrp(reducer[j],r);
527    }
528    poly hh=nc_rat_ReduceSpolyNew(reducer[j], *h, ishift, r);
529    //    p_Delete(h,r);
530    *h=hh;
531    if (TEST_OPT_DEBUG)
532    {
533      PrintS(" to ");
534      p_wrp(*h,r);
535      PrintLn();
536    }
537    if ((*h)==NULL)
538    {
539      return 0;
540    }
541  }
542}
543#endif
Note: See TracBrowser for help on using the repository browser.