source: git/kernel/kbuckets.cc @ 86016d

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