source: git/kernel/kbuckets.cc @ 811357f

spielwiese
Last change on this file since 811357f was 811357f, checked in by Michael Brickenstein <bricken@…>, 18 years ago
*bricken: prot git-svn-id: file:///usr/local/Singular/svn/trunk@8865 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: kbuckets.cc,v 1.11 2005-12-19 13:33:50 bricken Exp $ */
5
6#include "mod2.h"
7#include "structs.h"
8#include "omalloc.h"
9#include "p_polys.h"
10#include "febase.h"
11#include "pShallowCopyDelete.h"
12#include "kbuckets.h"
13#include "numbers.h"
14#include "p_Procs.h"
15
16#ifdef HAVE_COEF_BUCKETS
17#define USE_COEF_BUCKETS
18#endif
19
20#ifdef USE_COEF_BUCKETS
21#define MULTIPLY_BUCKET(B,I) do                                        \
22  { if (B->coef[I]!=NULL)                                              \
23    {                                                                  \
24      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
25      B->coef[I]=NULL;                                                 \
26    }                                                                  \
27  } while(0)
28#else
29#define MULTIPLY_BUCKET(B,I)
30#endif
31static omBin kBucket_bin = omGetSpecBin(sizeof(kBucket));
32static int coef_start=1;
33//////////////////////////////////////////////////////////////////////////
34///
35/// Some internal stuff
36///
37
38// returns ceil(log_4(l))
39inline unsigned int pLogLength(unsigned int l)
40{
41  unsigned int i = 0;
42
43  if (l == 0) return 0;
44  l--;
45#ifdef BUCKET_TWO_BASE
46  while ((l = (l >> 1))) i++;
47#else
48  while ((l = (l >> 2))) i++;
49#endif
50  return i+1;
51}
52
53// returns ceil(log_4(pLength(p)))
54inline unsigned int pLogLength(poly p)
55{
56  return pLogLength((unsigned int) pLength(p));
57}
58
59#ifdef KDEBUG
60
61#ifndef HAVE_PSEUDO_BUCKETS
62BOOLEAN kbTest_i(kBucket_pt bucket, int i)
63{//sBucketSortMerge
64  #ifdef USE_COEF_BUCKETS
65  assume(bucket->coef[0]==NULL);
66  if ((bucket->coef[i]!=NULL) && (bucket->buckets[i]==NULL))
67  {
68    dReportError("Bucket %d coef not NULL", i);
69  }
70  if (bucket->coef[i]!=NULL)
71  {
72    assume(bucket->buckets[i]!=NULL);
73    _p_Test(bucket->coef[i],bucket->bucket_ring,PDEBUG);
74    }
75  #endif
76  pFalseReturn(p_Test(bucket->buckets[i], bucket->bucket_ring));
77  if (bucket->buckets_length[i] != pLength(bucket->buckets[i]))
78  {
79    dReportError("Bucket %d lengths difference should:%d has:%d",
80                 i, bucket->buckets_length[i], pLength(bucket->buckets[i]));
81  }
82  else if (i > 0 && (int) pLogLength(bucket->buckets_length[i]) > i)
83  {
84    dReportError("Bucket %d too long %d",
85                 i, bucket->buckets_length[i]);
86  }
87  if (i==0 && bucket->buckets_length[0] > 1)
88  {
89    dReportError("Bucket 0 too long");
90  }
91  return TRUE;
92}
93
94
95BOOLEAN kbTest(kBucket_pt bucket)
96{
97  assume(bucket->coef[0]==NULL);
98  int i;
99  poly lm = bucket->buckets[0];
100
101  omCheckAddrBin(bucket, kBucket_bin);
102  if (! kbTest_i(bucket, 0)) return FALSE;
103  for (i=1; i<= (int) bucket->buckets_used; i++)
104  {
105    if (!kbTest_i(bucket, i)) return FALSE;
106    if (lm != NULL &&  bucket->buckets[i] != NULL
107        && p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) != 1)
108    {
109      dReportError("Bucket %d larger or equal than lm", i);
110      if (p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) ==0)
111      dReportError("Bucket %d equal to lm", i);
112      return FALSE;
113    }
114    if (!p_Test(bucket->buckets[i],bucket->bucket_ring))
115    {
116      dReportError("Bucket %d is not =0(4)", i);
117      return FALSE;
118    }
119  }
120
121  for (; i<=MAX_BUCKET; i++)
122  {
123    if (bucket->buckets[i] != NULL || bucket->buckets_length[i] != 0)
124    {
125      dReportError("Bucket %d not zero", i);
126      return FALSE;
127    }
128  }
129  for(i=0;i<=MAX_BUCKET;i++){
130    if (bucket->buckets[i]!=NULL){
131        int j;
132        for(j=i+1;j<=MAX_BUCKET;j++){
133            if (bucket->buckets[j]==bucket->buckets[i])
134                dReportError("Bucket %d %d equal", i,j);
135                return FALSE;
136        }
137    }
138    if (bucket->coef[i]!=NULL){
139        int j;
140        for(j=i+1;j<=MAX_BUCKET;j++){
141            if (bucket->coef[j]==bucket->coef[i])
142                dReportError("internal coef %d %d equal", i,j);
143                return FALSE;
144        }
145        for(j=0;j<=MAX_BUCKET;j++){
146            if (bucket->coef[j]==bucket->coef[i])
147                dReportError("internal coef %d equals bucket %d", i,j);
148                return FALSE;
149        }
150    }
151  }
152  return TRUE;
153}
154
155#else // HAVE_PSEUDO_BUCKETS
156BOOLEAN kbTest(kBucket_pt bucket)
157{
158  return TRUE;
159}
160#endif // ! HAVE_PSEUDO_BUCKETS
161#endif // KDEBUG
162
163//////////////////////////////////////////////////////////////////////////
164///
165/// Creation/Destruction of buckets
166///
167
168kBucket_pt kBucketCreate(ring bucket_ring)
169{
170  assume(bucket_ring != NULL);
171  kBucket_pt bucket = (kBucket_pt) omAlloc0Bin(kBucket_bin);
172  bucket->bucket_ring = bucket_ring;
173  return bucket;
174}
175void kBucketDestroy(kBucket_pt *bucket_pt)
176{
177  omFreeBin(*bucket_pt, kBucket_bin);
178  *bucket_pt = NULL;
179}
180
181
182void kBucketDeleteAndDestroy(kBucket_pt *bucket_pt)
183{
184  kBucket_pt bucket = *bucket_pt;
185  kbTest(bucket);
186  int i;
187  for (i=0; i<= bucket->buckets_used; i++)
188  {
189    if (bucket->buckets[i] != NULL)
190    {
191      p_Delete(&(bucket->buckets[i]), bucket->bucket_ring);
192#ifdef USE_COEF_BUCKETS
193      if (bucket->coef[i]!=NULL)
194        p_Delete(&(bucket->coef[i]), bucket->bucket_ring);
195#endif
196    }
197  }
198  omFreeBin(bucket, kBucket_bin);
199  *bucket_pt = NULL;
200}
201
202/////////////////////////////////////////////////////////////////////////////
203// Convertion from/to Bpolys
204//
205#ifndef HAVE_PSEUDO_BUCKETS
206
207inline void kBucketMergeLm(kBucket_pt bucket)
208{
209
210  if (bucket->buckets[0] != NULL)
211  {
212    poly lm = bucket->buckets[0];
213    int i = 1;
214#ifdef BUCKET_TWO_BASE
215    int l = 2;
216    while ( bucket->buckets_length[i] >= l)
217    {
218      i++;
219      l = l << 1;
220    }
221#else
222    int l = 4;
223    while ( bucket->buckets_length[i] >= l)
224    {
225      i++;
226      l = l << 2;
227    }
228#endif
229#if 0
230    MULTIPLY_BUCKET(bucket,i);
231    pNext(lm) = bucket->buckets[i];
232    bucket->buckets[i] = lm;
233    bucket->buckets_length[i]++;
234    assume(i <= bucket->buckets_used+1);
235    if (i > bucket->buckets_used)  bucket->buckets_used = i;
236    bucket->buckets[0] = NULL;
237    bucket->buckets_length[0] = 0;
238#endif
239    if (i > bucket->buckets_used)  bucket->buckets_used = i;
240    assume(i!=0);
241    if (bucket->buckets[i]!=NULL){
242       MULTIPLY_BUCKET(bucket,i);
243       pNext(lm) = bucket->buckets[i];
244       bucket->buckets[i] = lm;
245       bucket->buckets_length[i]++;
246       assume(i <= bucket->buckets_used+1);
247       
248    } else {
249      #if 1
250       assume(bucket->buckets[i]==NULL);
251       assume(bucket->coef[0]==NULL);
252       assume(pLength(lm)==1);
253       assume(pNext(lm)==NULL);
254       number coef=p_GetCoeff(lm,bucket->bucket_ring);
255       p_SetCoeff0(lm, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
256       bucket->buckets[i]=lm;
257       bucket->buckets_length[i]=1;
258       bucket->coef[i]=p_NSet(n_Copy(coef,bucket->bucket_ring),bucket->bucket_ring);
259       
260       bucket->buckets[i]=lm;
261       bucket->buckets_length[i]=1;
262       #else
263       MULTIPLY_BUCKET(bucket,i);
264       pNext(lm) = bucket->buckets[i];
265       bucket->buckets[i] = lm;
266       bucket->buckets_length[i]++;
267       assume(i <= bucket->buckets_used+1);
268       #endif
269    }
270    bucket->buckets[0]=NULL;
271    bucket->buckets_length[0] = 0;
272    bucket->coef[0]=NULL;
273    kbTest(bucket);
274  }
275}
276
277BOOLEAN kBucketIsCleared(kBucket_pt bucket)
278{
279  int i;
280
281  for (i = 0;i<=MAX_BUCKET;i++)
282  {
283    if (bucket->buckets[i] != NULL) return FALSE;
284    if (bucket->coef[i] != NULL) return FALSE;
285    if (bucket->buckets_length[i] != 0) return FALSE;
286  }
287  return TRUE;
288}
289
290void kBucketInit(kBucket_pt bucket, poly lm, int length)
291{
292  //assume(false);
293  assume(bucket != NULL);
294  assume(length <= 0 || length == pLength(lm));
295  assume(kBucketIsCleared(bucket));
296
297  if (lm == NULL) return;
298
299  if (length <= 0)
300    length = pLength(lm);
301
302  bucket->buckets[0] = lm;
303  assume(bucket->coef[0]==NULL);
304  bucket->coef[0]=NULL;
305  if (lm!=NULL)
306    bucket->buckets_length[0] = 1;
307  else
308    bucket->buckets_length[0]= 0;
309  if (length > 1)
310  {
311    unsigned int i = pLogLength(length-1);
312    bucket->buckets[i] = pNext(lm);
313    pNext(lm) = NULL;
314    bucket->buckets_length[i] = length-1;
315    bucket->buckets_used = i;
316  }
317  else
318  {
319    bucket->buckets_used = 0;
320  }
321}
322
323int kBucketCanonicalize(kBucket_pt bucket)
324{
325  MULTIPLY_BUCKET(bucket,1);
326  kbTest(bucket);
327  poly p = bucket->buckets[1];
328  poly lm;
329  int pl = bucket->buckets_length[1];//, i;
330  int i;
331  bucket->buckets[1] = NULL;
332  bucket->buckets_length[1] = 0;
333  assume(bucket->coef[1]==NULL);
334  ring r=bucket->bucket_ring;
335
336
337  for (i=1; i<=bucket->buckets_used; i++)
338  {
339  #ifdef USE_COEF_BUCKETS
340    if (bucket->coef[i]!=NULL)
341    {
342      assume(bucket->buckets[i]!=NULL);
343      p = p_Plus_mm_Mult_qq(p, bucket->coef[i], bucket->buckets[i],
344                 pl, bucket->buckets_length[i], r);
345      p_Delete(&bucket->coef[i],r);
346      p_Delete(&bucket->buckets[i],r);
347    }
348    else
349    p = p_Add_q(p, bucket->buckets[i],
350                 pl, bucket->buckets_length[i], r);
351  #else
352    p = p_Add_q(p, bucket->buckets[i],
353                 pl, bucket->buckets_length[i], r);
354  #endif
355    if (i==1) continue;
356    bucket->buckets[i] = NULL;
357    bucket->buckets_length[i] = 0;
358  }
359  assume(bucket->coef[0]==NULL);
360  lm = bucket->buckets[0];
361  if (lm != NULL)
362  {
363    pNext(lm) = p;
364    p = lm;
365    pl++;
366    bucket->buckets[0] = NULL;
367    bucket->buckets_length[0] = 0;
368  }
369  if (pl > 0)
370  {
371    i = pLogLength(pl);
372    bucket->buckets[i] = p;
373    bucket->buckets_length[i] = pl;
374  }
375  else
376  {
377    i = 0;
378  }
379  bucket->buckets_used = i;
380  assume(bucket->coef[0]==NULL);
381  assume(bucket->coef[i]==NULL);
382  assume(pLength(p) == (int) pl);
383  kbTest(bucket);
384  return i;
385}
386
387void kBucketClear(kBucket_pt bucket, poly *p, int *length)
388{
389  int i = kBucketCanonicalize(bucket);
390  if (i > 0)
391  {
392  #ifdef USE_COEF_BUCKETS
393    MULTIPLY_BUCKET(bucket,i);
394    //bucket->coef[i]=NULL;
395#endif
396    *p = bucket->buckets[i];
397    *length = bucket->buckets_length[i];
398    bucket->buckets[i] = NULL;
399    bucket->buckets_length[i] = 0;
400    bucket->buckets_used = 0;
401
402  }
403  else
404  {
405    *p = NULL;
406    *length = 0;
407  }
408}
409
410void kBucketSetLm(kBucket_pt bucket, poly lm)
411{
412  kBucketMergeLm(bucket);
413  pNext(lm) = NULL;
414  bucket->buckets[0] = lm;
415  bucket->buckets_length[0] = 1;
416}
417
418#else // HAVE_PSEUDO_BUCKETS
419
420void kBucketInit(kBucket_pt bucket, poly lm, int length)
421{
422  int i;
423
424  assume(bucket != NULL);
425  assume(length <= 0 || length == pLength(lm));
426
427  bucket->p = lm;
428  if (length <= 0) bucket->l = pLength(lm);
429  else bucket->l = length;
430
431}
432
433const poly kBucketGetLm(kBucket_pt bucket)
434{
435  return bucket->p;
436}
437
438poly kBucketExtractLm(kBucket_pt bucket)
439{
440  poly lm = bucket->p;
441  assume(pLength(bucket->p) == bucket->l);
442  pIter(bucket->p);
443  (bucket->l)--;
444  pNext(lm) = NULL;
445  return lm;
446}
447
448void kBucketClear(kBucket_pt bucket, poly *p, int *length)
449{
450  assume(pLength(bucket->p) == bucket->l);
451  *p = bucket->p;
452  *length = bucket->l;
453  bucket->p = NULL;
454  bucket->l = 0;
455}
456
457#endif // ! HAVE_PSEUDO_BUCKETS
458//////////////////////////////////////////////////////////////////////////
459///
460/// For changing the ring of the Bpoly to new_tailBin
461///
462void kBucketShallowCopyDelete(kBucket_pt bucket,
463                              ring new_tailRing, omBin new_tailBin,
464                              pShallowCopyDeleteProc p_shallow_copy_delete)
465{
466#ifndef HAVE_PSEUDO_BUCKETS
467  int i;
468
469  kBucketCanonicalize(bucket);
470  for (i=0; i<= bucket->buckets_used; i++)
471    if (bucket->buckets[i] != NULL)
472    {
473      MULTIPLY_BUCKET(bucket,i);
474      bucket->buckets[i] = p_shallow_copy_delete(bucket->buckets[i],
475                                                 bucket->bucket_ring,
476                                                 new_tailRing,
477                                                 new_tailBin);
478    }
479#else
480  bucket->p = p_shallow_copy_delete(p,
481                                    bucket_ring,
482                                    new_tailRing,
483                                    new_tailBin);
484#endif
485  bucket->bucket_ring = new_tailRing;
486}
487
488
489
490//////////////////////////////////////////////////////////////////////////
491///
492/// Multiply Bucket by number ,i.e. Bpoly == n*Bpoly
493///
494void kBucket_Mult_n(kBucket_pt bucket, number n)
495{
496#ifndef HAVE_PSEUDO_BUCKETS
497  kbTest(bucket);
498  ring r=bucket->bucket_ring;
499  int i;
500
501  for (i=0; i<= bucket->buckets_used; i++)
502  {
503    if (bucket->buckets[i] != NULL)
504    {
505#ifdef USE_COEF_BUCKETS
506      if (i<coef_start)
507        bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
508      else
509      if (bucket->coef[i]!=NULL)
510      {
511        bucket->coef[i] = p_Mult_nn(bucket->coef[i],n,r);
512      }
513      else
514      {
515        bucket->coef[i] = p_NSet(n_Copy(n,r),r);
516      }
517#else
518      bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
519#endif
520    }
521  }
522  kbTest(bucket);
523#else
524  bucket->p = p_Mult_nn(bucket->p, n, bucket->bucket_ring);
525#endif
526}
527
528
529//////////////////////////////////////////////////////////////////////////
530///
531/// Add to Bucket a poly ,i.e. Bpoly == q+Bpoly
532///
533void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
534{
535  if (q == NULL) return;
536  assume(*l <= 0 || pLength(q) == *l);
537
538  int i, l1;
539  ring r = bucket->bucket_ring;
540
541  if (*l <= 0)
542  {
543    l1 = pLength(q);
544    *l = l1;
545  }
546  else
547    l1 = *l;
548
549  kBucketMergeLm(bucket);
550  kbTest(bucket);
551  i = pLogLength(l1);
552
553  while (bucket->buckets[i] != NULL)
554  {
555    //MULTIPLY_BUCKET(bucket,i);
556  #ifdef USE_COEF_BUCKETS
557    if (bucket->coef[i]!=NULL)
558    {
559      q = p_Plus_mm_Mult_qq(q, bucket->coef[i], bucket->buckets[i],
560                 l1, bucket->buckets_length[i], r);
561      p_Delete(&bucket->coef[i],r);
562      p_Delete(&bucket->buckets[i],r);
563    }
564    else
565    q = p_Add_q(q, bucket->buckets[i],
566                 l1, bucket->buckets_length[i], r);
567  #else
568    q = p_Add_q(q, bucket->buckets[i],
569                 l1, bucket->buckets_length[i], r);
570  #endif
571    bucket->buckets[i] = NULL;
572    bucket->buckets_length[i] = 0;
573    i = pLogLength(l1);
574  }
575
576  bucket->buckets[i] = q;
577  bucket->buckets_length[i]=l1;
578  if (i >= bucket->buckets_used)
579    bucket->buckets_used = i;
580  else
581    kBucketAdjustBucketsUsed(bucket);
582  kbTest(bucket);
583}
584
585
586
587//////////////////////////////////////////////////////////////////////////
588///
589/// Bpoly == Bpoly - m*p; where m is a monom
590/// Does not destroy p and m
591/// assume (*l <= 0 || pLength(p) == *l)
592void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l,
593                            poly spNoether)
594{
595  assume(*l <= 0 || pLength(p) == *l);
596  int i, l1;
597  poly p1 = p;
598  poly last;
599  ring r = bucket->bucket_ring;
600
601  if (*l <= 0)
602  {
603    l1 = pLength(p1);
604    *l = l1;
605  }
606  else
607    l1 = *l;
608
609  if (m == NULL || p == NULL) return;
610
611#ifndef HAVE_PSEUDO_BUCKETS
612  kBucketMergeLm(bucket);
613  kbTest(bucket);
614  i = pLogLength(l1);
615
616  if ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
617  {
618    assume(pLength(bucket->buckets[i])==bucket->buckets_length[i]);
619//#ifdef USE_COEF_BUCKETS
620//     if(bucket->coef[i]!=NULL)
621//     {
622//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
623//       bucket->coef[i]=NULL;
624//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
625//                               bucket->buckets_length[i], l1,
626//                             spNoether, r);
627//     }
628//     else
629//#endif
630    MULTIPLY_BUCKET(bucket,i);
631    p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
632                            bucket->buckets_length[i], l1,
633                            spNoether, r);
634    l1 = bucket->buckets_length[i];
635    bucket->buckets[i] = NULL;
636    bucket->buckets_length[i] = 0;
637    i = pLogLength(l1);
638  }
639  else
640  {
641    pSetCoeff0(m, nNeg(pGetCoeff(m)));
642    if (spNoether != NULL)
643    {
644      l1 = -1;
645      p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r, last);
646      i = pLogLength(l1);
647    }
648    else
649      p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
650    pSetCoeff0(m, nNeg(pGetCoeff(m)));
651  }
652
653  while (bucket->buckets[i] != NULL)
654  {
655    //kbTest(bucket);
656    MULTIPLY_BUCKET(bucket,i);
657    p1 = p_Add_q(p1, bucket->buckets[i],
658                 l1, bucket->buckets_length[i], r);
659    bucket->buckets[i] = NULL;
660    bucket->buckets_length[i] = 0;
661    i = pLogLength(l1);
662  }
663
664  bucket->buckets[i] = p1;
665  bucket->buckets_length[i]=l1;
666  if (i >= bucket->buckets_used)
667    bucket->buckets_used = i;
668  else
669    kBucketAdjustBucketsUsed(bucket);
670#else // HAVE_PSEUDO_BUCKETS
671  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
672                               bucket->l, l1,
673                               spNoether, r);
674#endif
675}
676
677//////////////////////////////////////////////////////////////////////////
678///
679/// Bpoly == Bpoly + m*p; where m is a monom
680/// Does not destroy p and m
681/// assume (l <= 0 || pLength(p) == l)
682void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
683{
684  assume(l <= 0 || pLength(p) == l);
685  int i, l1;
686  poly p1 = p;
687  poly last;
688  ring r = bucket->bucket_ring;
689
690  if (l <= 0)
691  {
692    l1 = pLength(p1);
693    l = l1;
694  }
695  else
696    l1 = l;
697
698  if (m == NULL || p == NULL) return;
699
700  kBucketMergeLm(bucket);
701  kbTest(bucket);
702  i = pLogLength(l1);
703  number n=n_Init(1,r);
704  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
705  { 
706    //if (FALSE){
707    if ((bucket->coef[i]!=NULL) &&(i>=coef_start)){
708        number orig_coef=p_GetCoeff(bucket->coef[i],r);
709        //we take ownership:
710        p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
711        number add_coef=n_Copy(p_GetCoeff(m,r),r);
712        number gcd=n_Gcd(add_coef, orig_coef,r);
713       
714        if (!(n_IsOne(gcd,r))) {
715            number orig_coef2=n_IntDiv(orig_coef,gcd,r);
716            number add_coef2=n_IntDiv(add_coef, gcd,r);
717            n_Delete(&orig_coef,r);
718            n_Delete(&add_coef,r);
719            orig_coef=orig_coef2;
720            add_coef=add_coef2;
721         
722            //p_Mult_nn(bucket->buckets[i], orig_coef,r);
723            n_Delete(&n,r);
724            n=gcd;
725           
726        }
727       
728        //assume(n_IsOne(n,r));
729        number backup=p_GetCoeff(m,r);
730       
731        p_SetCoeff0(m,add_coef,r);
732        bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
733       
734        n_Delete(&orig_coef,r);
735        p_Delete(&bucket->coef[i],r);
736       
737
738        p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
739                           bucket->buckets_length[i], l1, r);
740        bucket->buckets[i]=NULL;
741        l1=bucket->buckets_length[i];
742        assume(l1==pLength(p1));
743       
744        p_SetCoeff(m,backup,r); //deletes add_coef
745       
746       
747    }
748    else {
749   
750    MULTIPLY_BUCKET(bucket,i);
751    p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
752                           bucket->buckets_length[i], l1, r);
753    l1 = bucket->buckets_length[i];
754    bucket->buckets[i] = NULL;
755    bucket->buckets_length[i] = 0;
756    i = pLogLength(l1);
757    }
758  }
759  else
760  {
761   
762    number swap_n=p_GetCoeff(m,r);
763   
764    assume(n_IsOne(n,r));
765    p_SetCoeff0(m,n,r);
766    n=swap_n;
767    //p_SetCoeff0(n, swap_n, r);
768    //p_GetCoeff0(n, swap_n,r);
769   
770    p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
771    //m may not be changed
772    p_SetCoeff(m,n_Copy(n,r),r);
773  }
774
775  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
776  {
777
778    assume(i!=0);
779    if ((bucket->coef[i]!=NULL) &&(i>=coef_start)){
780        number orig_coef=p_GetCoeff(bucket->coef[i],r);
781        //we take ownership:
782        p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
783        number add_coef=n;
784        number gcd=n_Gcd(add_coef, orig_coef,r);
785       
786        if (!(n_IsOne(gcd,r))) {
787            number orig_coef2=n_IntDiv(orig_coef,gcd,r);
788            number add_coef2=n_IntDiv(add_coef, gcd,r);
789            n_Delete(&orig_coef,r);
790            n_Delete(&add_coef,r);
791            orig_coef=orig_coef2;
792            add_coef=add_coef2;
793            //p_Mult_nn(bucket->buckets[i], orig_coef,r);
794            n_Delete(&n,r);
795            n=gcd;
796        }
797       
798        //assume(n_IsOne(n,r));
799        bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
800        p1=p_Mult_nn(p1,add_coef,r);
801       
802     
803       
804
805        p1 = p_Add_q(p1, bucket->buckets[i],r);
806        l1=pLength(p1);
807       
808        bucket->buckets[i]=NULL;
809        n_Delete(&orig_coef,r);
810        p_Delete(&bucket->coef[i],r);
811        //l1=bucket->buckets_length[i];
812        assume(l1==pLength(p1));
813       
814
815       
816    } else {   
817   
818    //don't do that, pull out gcd
819    if(!(n_IsOne(n,r))) {
820        p1=p_Mult_nn(p1, n, r);
821        n_Delete(&n,r);
822        n=n_Init(1,r);
823    } 
824   
825    MULTIPLY_BUCKET(bucket,i);
826    p1 = p_Add_q(p1, bucket->buckets[i],
827                 l1, bucket->buckets_length[i], r);
828    bucket->buckets[i] = NULL;
829    bucket->buckets_length[i] = 0;
830    i = pLogLength(l1);
831    }
832  }
833
834 
835  bucket->buckets[i] = p1;
836 
837  assume(bucket->coef[i]==NULL);
838  if (!(n_IsOne(n,r))){
839    bucket->coef[i]=p_NSet(n,r);
840  } else {
841    bucket->coef[i]=NULL;
842    n_Delete(&n,r);
843  }
844  if ((p1==NULL) && (bucket->coef[i]!=NULL)) p_Delete(&bucket->coef[i],r);
845  bucket->buckets_length[i]=l1;
846  if (i >= bucket->buckets_used)
847    bucket->buckets_used = i;
848  else
849    kBucketAdjustBucketsUsed(bucket);
850
851  kbTest(bucket);
852}
853
854poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
855{
856  if (q == NULL) return append;
857  poly lm;
858  loop
859  {
860    lm = kBucketGetLm(bucket);
861    if (lm == NULL) return append;
862    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
863    {
864      lm = kBucketExtractLm(bucket);
865      pNext(append) = lm;
866      pIter(append);
867    }
868    else
869    {
870      return append;
871    }
872  }
873}
874
875/////////////////////////////////////////////////////////////////////////////
876//
877// Extract all monomials from bucket with component comp
878// Return as a polynomial *p with length *l
879// In other words, afterwards
880// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
881// and components of monomials of *p are all 0
882//
883
884// Hmm... for now I'm too lazy to implement those independent of currRing
885// But better declare it extern than including polys.h
886extern void pTakeOutComp(poly *p, Exponent_t comp, poly *q, int *lq);
887void pDecrOrdTakeOutComp(poly *p, Exponent_t comp, Order_t order,
888                         poly *q, int *lq);
889
890void kBucketTakeOutComp(kBucket_pt bucket,
891                        Exponent_t comp,
892                        poly *r_p, int *l)
893{
894  poly p = NULL, q;
895  int i, lp = 0, lq;
896
897#ifndef HAVE_PSEUDO_BUCKETS
898  kBucketMergeLm(bucket);
899  for (i=1; i<=bucket->buckets_used; i++)
900  {
901    if (bucket->buckets[i] != NULL)
902    {
903      MULTIPLY_BUCKET(bucket,i);
904      pTakeOutComp(&(bucket->buckets[i]), comp, &q, &lq);
905      if (q != NULL)
906      {
907        assume(pLength(q) == lq);
908        bucket->buckets_length[i] -= lq;
909        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
910        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
911      }
912    }
913  }
914  kBucketAdjustBucketsUsed(bucket);
915#else
916  pTakeOutComp(&(bucket->p), comp, &p, &lp);
917  (bucket->l) -= lp;
918#endif
919  *r_p = p;
920  *l = lp;
921
922  kbTest(bucket);
923}
924
925void kBucketDecrOrdTakeOutComp(kBucket_pt bucket,
926                               Exponent_t comp, Order_t order,
927                               poly *r_p, int *l)
928{
929  poly p = NULL, q;
930  int i, lp = 0, lq;
931
932#ifndef HAVE_PSEUDO_BUCKETS
933  kBucketMergeLm(bucket);
934  for (i=1; i<=bucket->buckets_used; i++)
935  {
936    if (bucket->buckets[i] != NULL)
937    {
938      MULTIPLY_BUCKET(bucket,i);
939      pDecrOrdTakeOutComp(&(bucket->buckets[i]), comp, order, &q, &lq);
940      if (q != NULL)
941      {
942        bucket->buckets_length[i] -= lq;
943        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
944      }
945    }
946  }
947  kBucketAdjustBucketsUsed(bucket);
948#else
949  pDecrOrdTakeOutComp(&(bucket->p), comp, order, &p, &lp);
950  (bucket->l) -= lp;
951#endif
952
953  *r_p = p;
954  *l = lp;
955}
956
957/////////////////////////////////////////////////////////////////////////////
958// Reduction of Bpoly with a given poly
959//
960
961extern int ksCheckCoeff(number *a, number *b);
962
963number kBucketPolyRed(kBucket_pt bucket,
964                      poly p1, int l1,
965                      poly spNoether)
966{
967  assume(p1 != NULL &&
968         p_DivisibleBy(p1,  kBucketGetLm(bucket), bucket->bucket_ring));
969  assume(pLength(p1) == (int) l1);
970
971  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
972  BOOLEAN reset_vec=FALSE;
973  number rn;
974
975  if(a1==NULL)
976  {
977    p_DeleteLm(&lm, bucket->bucket_ring);
978    return nInit(1);
979  }
980
981  if (! nIsOne(pGetCoeff(p1)))
982  {
983    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
984    int ct = ksCheckCoeff(&an, &bn);
985    p_SetCoeff(lm, bn, bucket->bucket_ring);
986    if ((ct == 0) || (ct == 2)) kBucket_Mult_n(bucket, an);
987    rn = an;
988  }
989  else
990  {
991    rn = nInit(1);
992  }
993
994  if (p_GetComp(p1, bucket->bucket_ring) != p_GetComp(lm, bucket->bucket_ring))
995  {
996    p_SetCompP(a1, p_GetComp(lm, bucket->bucket_ring), bucket->bucket_ring);
997    reset_vec = TRUE;
998    p_SetComp(lm, p_GetComp(p1, bucket->bucket_ring), bucket->bucket_ring);
999    p_Setm(lm, bucket->bucket_ring);
1000  }
1001
1002  p_ExpVectorSub(lm,p1, bucket->bucket_ring);
1003  l1--;
1004
1005  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1006
1007  p_DeleteLm(&lm, bucket->bucket_ring);
1008  if (reset_vec) p_SetCompP(a1, 0, bucket->bucket_ring);
1009  kbTest(bucket);
1010  return rn;
1011}
1012static BOOLEAN nIsPseudoUnit(number n, ring r){
1013    if (rField_is_Zp(r))
1014        return TRUE;
1015       
1016    //if (r->parameter!=NULL)
1017    number one=n_Init(1,r);
1018    if (n_Equal(n,one,r)) {
1019    n_Delete(&one,r);
1020    return TRUE;
1021    }
1022    n_Delete(&one,r);
1023    number minus_one=n_Init(-1,r);
1024    if (n_Equal(n,minus_one,r)){
1025        n_Delete(&minus_one,r);
1026        return TRUE;
1027    }
1028
1029    return FALSE;
1030   
1031}
1032void kBucketSimpleContent(kBucket_pt bucket){
1033    int i;
1034    //PrintS("HHHHHHHHHHHHH");
1035    for (i=0;i<=MAX_BUCKET;i++){
1036        //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1037        //    PrintS("H2H2H2");
1038          if (i==0)
1039           
1040           {
1041           
1042           assume(bucket->buckets[i]==NULL);
1043           }
1044        if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL)){
1045            return;
1046         }
1047    }
1048    //return;
1049    ring r=bucket->bucket_ring;
1050    number coef=n_Init(0,r);
1051    //ATTENTION: will not work correct for GB over ring
1052    //if (TEST_OPT_PROT)
1053    //    PrintS("CCCCCCCCCCCCC");
1054    for (i=0;i<=MAX_BUCKET;i++){
1055        if (i==0)
1056           
1057           {
1058           
1059           assume(bucket->buckets[i]==NULL);
1060           }
1061        if (bucket->buckets[i]!=NULL){
1062           
1063            assume(bucket->coef[i]!=NULL);
1064            //in this way it should crash on programming errors, yeah
1065            number temp=nGcd(coef, pGetCoeff(bucket->buckets[i]),r);
1066            n_Delete(&coef,r );
1067            coef=temp;
1068            if (nIsPseudoUnit(coef,r))
1069            {
1070                n_Delete(&coef,r);
1071                return;
1072            }
1073           
1074         }
1075    }
1076    if (TEST_OPT_PROT)
1077        PrintS("S");
1078    for(i=0;i<=MAX_BUCKET;i++){
1079        if (bucket->buckets[i]!=NULL)
1080        {
1081        assume(bucket->coef[i]!=NULL);
1082        number lc=p_GetCoeff(bucket->coef[i],r);
1083        p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1084        }
1085    }
1086    return;
1087}
1088
1089
1090poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i){
1091    assume(bucket->buckets[i]!=NULL);
1092   
1093    ring r=bucket->bucket_ring;
1094    poly p=bucket->buckets[i];
1095    bucket->buckets_length[i]--;
1096    if (bucket->coef[i]!=NULL){
1097        poly next=pNext(p);
1098        if (next==NULL){
1099            MULTIPLY_BUCKET(bucket,i);
1100            p=bucket->buckets[i];
1101            bucket->buckets[i]=NULL;
1102            return p;
1103        } else {
1104            bucket->buckets[i]=next;
1105            number c=p_GetCoeff(bucket->coef[i],r);
1106            pNext(p)=NULL;
1107            p=p_Mult_nn(p,c,r);
1108            assume(p!=NULL);
1109            return p;
1110        }
1111       
1112    } else {
1113       
1114        bucket->buckets[i]=pNext(bucket->buckets[i]);
1115        pNext(p)=NULL;
1116        assume(p!=NULL);
1117        return p;
1118    }
1119   
1120}
Note: See TracBrowser for help on using the repository browser.