source: git/kernel/kbuckets.cc @ 4ecd42

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