source: git/libpolys/polys/kbuckets.cc @ d5a330

jengelh-datetimespielwiese
Last change on this file since d5a330 was d5a330, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: unify changes for HAVE_RING and HAVE_PLURAL
  • Property mode set to 100644
File size: 32.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5//#include <kernel/mod2.h>
6#include <omalloc/omalloc.h>
7#include <polys/monomials/p_polys.h>
8//#include <kernel/febase.h>
9//#include <kernel/pShallowCopyDelete.h>
10#include <coeffs/coeffs.h>
11#include <polys/monomials/ring.h>
12//#include <kernel/p_Procs.h>
13//#include <kernel/kutil.h>
14#include <polys/kbuckets.h>
15
16// #include <polys/operations/pShallowCopyDelete.h>
17
18
19#ifdef HAVE_COEF_BUCKETS
20#define USE_COEF_BUCKETS
21#endif
22
23#ifdef USE_COEF_BUCKETS
24#ifdef HAVE_RINGS_OLD
25#define MULTIPLY_BUCKET(B,I) do                                        \
26  { if (B->coef[I]!=NULL)                                              \
27    {                                                                  \
28      assume(p_IsConstant(B->Coef[i],B->bucket->ring));           \
29      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
30      B->coef[I]=NULL;                                                 \
31    }                                                                  \
32  } while(0)                                                           \
33    if (rField_is_Ring(B->bucket_ring)) B->buckets_length[i] = pLength(B->buckets[i]);
34#else
35#define MULTIPLY_BUCKET(B,I) do                                        \
36  { if (B->coef[I]!=NULL)                                              \
37    {                                                                  \
38      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
39      B->coef[I]=NULL;                                                 \
40    }                                                                  \
41  } while(0)
42#endif
43#else
44#define MULTIPLY_BUCKET(B,I)
45#endif
46static omBin kBucket_bin = omGetSpecBin(sizeof(kBucket));
47#ifdef USE_COEF_BUCKETS
48static int coef_start=1;
49#endif
50//////////////////////////////////////////////////////////////////////////
51///
52/// Some internal stuff
53///
54
55// returns ceil(log_4(l))
56inline unsigned int pLogLength(unsigned int l)
57{
58  unsigned int i = 0;
59
60  if (l == 0) return 0;
61  l--;
62#ifdef BUCKET_TWO_BASE
63  while ((l = (l >> 1))) i++;
64#else
65  while ((l = (l >> 2))) i++;
66#endif
67  return i+1;
68}
69
70// returns ceil(log_4(pLength(p)))
71inline unsigned int pLogLength(poly p)
72{
73  return pLogLength((unsigned int) pLength(p));
74}
75
76#ifdef KDEBUG
77
78#ifndef HAVE_PSEUDO_BUCKETS
79BOOLEAN kbTest_i(kBucket_pt bucket, int i)
80{//sBucketSortMerge
81  #ifdef USE_COEF_BUCKETS
82  assume(bucket->coef[0]==NULL);
83  if ((bucket->coef[i]!=NULL) && (bucket->buckets[i]==NULL))
84  {
85    dReportError("Bucket %d coef not NULL", i);
86  }
87  if (bucket->coef[i]!=NULL)
88  {
89    assume(bucket->buckets[i]!=NULL);
90    _p_Test(bucket->coef[i],bucket->bucket_ring,PDEBUG);
91    }
92  #endif
93  pFalseReturn(p_Test(bucket->buckets[i], bucket->bucket_ring));
94  if (bucket->buckets_length[i] != pLength(bucket->buckets[i]))
95  {
96    dReportError("Bucket %d lengths difference should:%d has:%d",
97                 i, bucket->buckets_length[i], pLength(bucket->buckets[i]));
98  }
99  else if (i > 0 && (int) pLogLength(bucket->buckets_length[i]) > i)
100  {
101    dReportError("Bucket %d too long %d",
102                 i, bucket->buckets_length[i]);
103  }
104  if (i==0 && bucket->buckets_length[0] > 1)
105  {
106    dReportError("Bucket 0 too long");
107  }
108  return TRUE;
109}
110
111
112BOOLEAN kbTest(kBucket_pt bucket)
113{
114  #ifdef HAVE_COEF_BUCKETS
115  assume(bucket->coef[0]==NULL);
116  #endif
117  int i;
118  poly lm = bucket->buckets[0];
119
120  omCheckAddrBin(bucket, kBucket_bin);
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  //if (TEST_OPT_PROT) { Print("C(%d)",pl); }
431  kbTest(bucket);
432  return i;
433}
434
435void kBucketClear(kBucket_pt bucket, poly *p, int *length)
436{
437  int i = kBucketCanonicalize(bucket);
438  if (i > 0)
439  {
440  #ifdef USE_COEF_BUCKETS
441    MULTIPLY_BUCKET(bucket,i);
442    //bucket->coef[i]=NULL;
443  #endif
444    *p = bucket->buckets[i];
445    *length = bucket->buckets_length[i];
446    bucket->buckets[i] = NULL;
447    bucket->buckets_length[i] = 0;
448    bucket->buckets_used = 0;
449
450  }
451  else
452  {
453    *p = NULL;
454    *length = 0;
455  }
456}
457
458void kBucketSetLm(kBucket_pt bucket, poly lm)
459{
460  kBucketMergeLm(bucket);
461  pNext(lm) = NULL;
462  bucket->buckets[0] = lm;
463  bucket->buckets_length[0] = 1;
464}
465
466#else // HAVE_PSEUDO_BUCKETS
467
468void kBucketInit(kBucket_pt bucket, poly lm, int length)
469{
470  int i;
471
472  assume(bucket != NULL);
473  assume(length <= 0 || length == pLength(lm));
474
475  bucket->p = lm;
476  if (length <= 0) bucket->l = pLength(lm);
477  else bucket->l = length;
478
479}
480
481const poly kBucketGetLm(kBucket_pt bucket)
482{
483  return bucket->p;
484}
485
486poly kBucketExtractLm(kBucket_pt bucket)
487{
488  poly lm = bucket->p;
489  assume(pLength(bucket->p) == bucket->l);
490  pIter(bucket->p);
491  (bucket->l)--;
492  pNext(lm) = NULL;
493  return lm;
494}
495
496void kBucketClear(kBucket_pt bucket, poly *p, int *length)
497{
498  assume(pLength(bucket->p) == bucket->l);
499  *p = bucket->p;
500  *length = bucket->l;
501  bucket->p = NULL;
502  bucket->l = 0;
503}
504
505#endif // ! HAVE_PSEUDO_BUCKETS
506//////////////////////////////////////////////////////////////////////////
507///
508/// For changing the ring of the Bpoly to new_tailBin
509///
510void kBucketShallowCopyDelete(kBucket_pt bucket,
511                              ring new_tailRing, omBin new_tailBin,
512                              pShallowCopyDeleteProc p_shallow_copy_delete)
513{
514#ifndef HAVE_PSEUDO_BUCKETS
515  int i;
516
517  kBucketCanonicalize(bucket);
518  for (i=0; i<= bucket->buckets_used; i++)
519    if (bucket->buckets[i] != NULL)
520    {
521      MULTIPLY_BUCKET(bucket,i);
522      bucket->buckets[i] = p_shallow_copy_delete(bucket->buckets[i],
523                                                 bucket->bucket_ring,
524                                                 new_tailRing,
525                                                 new_tailBin);
526    }
527#else
528  bucket->p = p_shallow_copy_delete(p,
529                                    bucket_ring,
530                                    new_tailRing,
531                                    new_tailBin);
532#endif
533  bucket->bucket_ring = new_tailRing;
534}
535
536//////////////////////////////////////////////////////////////////////////
537///
538/// Bucket number i from bucket is out of length sync, resync
539///
540void kBucketAdjust(kBucket_pt bucket, int i) {
541
542  MULTIPLY_BUCKET(bucket,i);
543
544  int l1 = bucket->buckets_length[i];
545  poly p1 = bucket->buckets[i];
546  bucket->buckets[i] = NULL;
547  bucket->buckets_length[i] = 0;
548  i = pLogLength(l1);
549
550  while (bucket->buckets[i] != NULL)
551  {
552    //kbTest(bucket);
553    MULTIPLY_BUCKET(bucket,i);
554    p1 = p_Add_q(p1, bucket->buckets[i],
555                 l1, bucket->buckets_length[i], bucket->bucket_ring);
556    bucket->buckets[i] = NULL;
557    bucket->buckets_length[i] = 0;
558    i = pLogLength(l1);
559  }
560
561  bucket->buckets[i] = p1;
562  bucket->buckets_length[i]=l1;
563  if (i >= bucket->buckets_used)
564    bucket->buckets_used = i;
565  else
566    kBucketAdjustBucketsUsed(bucket);
567}
568
569//////////////////////////////////////////////////////////////////////////
570///
571/// Multiply Bucket by number ,i.e. Bpoly == n*Bpoly
572///
573void kBucket_Mult_n(kBucket_pt bucket, number n)
574{
575#ifndef HAVE_PSEUDO_BUCKETS
576  kbTest(bucket);
577  ring r=bucket->bucket_ring;
578  int i;
579
580  for (i=0; i<= bucket->buckets_used; i++)
581  {
582    if (bucket->buckets[i] != NULL)
583    {
584#ifdef USE_COEF_BUCKETS
585      if (i<coef_start)
586        bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
587#ifdef HAVE_RINGS
588        /* Frank Seelisch on March 11, 2010:
589           This looks a bit strange: The following "if" is indented
590           like the previous line of code. But coded as it is,
591           it should actually be two spaces less indented.
592           Question: Should the following "if" also only be
593           performed when "(i<coef_start)" is true?
594           For the time being, I leave it as it is. */
595        if (rField_is_Ring(r) && !(rField_is_Domain(r)))
596        {
597          bucket->buckets_length[i] = pLength(bucket->buckets[i]);
598          kBucketAdjust(bucket, i);
599        }
600#endif
601      else
602      if (bucket->coef[i]!=NULL)
603      {
604        bucket->coef[i] = p_Mult_nn(bucket->coef[i],n,r);
605      }
606      else
607      {
608        bucket->coef[i] = p_NSet(n_Copy(n,r),r);
609      }
610#else
611      bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
612#ifdef HAVE_RINGS
613      if (rField_is_Ring(r) && !(rField_is_Domain(r)))
614      {
615        bucket->buckets_length[i] = pLength(bucket->buckets[i]);
616        kBucketAdjust(bucket, i);
617      }
618#endif
619#endif
620    }
621  }
622  kbTest(bucket);
623#else
624  bucket->p = p_Mult_nn(bucket->p, n, bucket->bucket_ring);
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  ring r = bucket->bucket_ring;
702
703  if (*l <= 0)
704  {
705    l1 = pLength(p1);
706    *l = l1;
707  }
708  else
709    l1 = *l;
710
711  if (m == NULL || p == NULL) return;
712
713#ifndef HAVE_PSEUDO_BUCKETS
714  kBucketMergeLm(bucket);
715  kbTest(bucket);
716  i = pLogLength(l1);
717
718#if defined(HAVE_RINGS)||defined(HAVE_PLURAL)
719  if ((rField_is_Ring(r) && !(rField_is_Domain(r)))
720  ||(rIsPluralRing(r)))
721  {
722    pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
723    p1=pp_Mult_mm(p,m,r);
724    pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
725    l1=pLength(p1);
726    i = pLogLength(l1);
727  }
728  else
729#endif
730  {
731    if ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
732    {
733      assume(pLength(bucket->buckets[i])==bucket->buckets_length[i]);
734//#ifdef USE_COEF_BUCKETS
735//     if(bucket->coef[i]!=NULL)
736//     {
737//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
738//       bucket->coef[i]=NULL;
739//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
740//                               bucket->buckets_length[i], l1,
741//                             spNoether, r);
742//     }
743//     else
744//#endif
745      MULTIPLY_BUCKET(bucket,i);
746      p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
747                            bucket->buckets_length[i], l1,
748                            spNoether, r);
749      l1 = bucket->buckets_length[i];
750      bucket->buckets[i] = NULL;
751      bucket->buckets_length[i] = 0;
752      i = pLogLength(l1);
753    }
754    else
755    {
756      pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
757      if (spNoether != NULL)
758      {
759        l1 = -1;
760        p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r);
761        i = pLogLength(l1);
762      }
763      else
764      {
765        p1 = r->p_Procs->pp_Mult_mm(p1, m, r);
766      }
767      pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
768    }
769  }
770
771  while (bucket->buckets[i] != NULL)
772  {
773    //kbTest(bucket);
774    MULTIPLY_BUCKET(bucket,i);
775    p1 = p_Add_q(p1, bucket->buckets[i],
776                 l1, bucket->buckets_length[i], r);
777    bucket->buckets[i] = NULL;
778    bucket->buckets_length[i] = 0;
779    i = pLogLength(l1);
780  }
781
782  bucket->buckets[i] = p1;
783  bucket->buckets_length[i]=l1;
784  if (i >= bucket->buckets_used)
785    bucket->buckets_used = i;
786  else
787    kBucketAdjustBucketsUsed(bucket);
788#else // HAVE_PSEUDO_BUCKETS
789  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
790                               bucket->l, l1,
791                               spNoether, r);
792#endif
793}
794
795//////////////////////////////////////////////////////////////////////////
796///
797/// Bpoly == Bpoly + m*p; where m is a monom
798/// Does not destroy p and m
799/// assume (l <= 0 || pLength(p) == l)
800void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
801{
802    assume((!rIsPluralRing(bucket->bucket_ring))||p_IsConstant(m, bucket->bucket_ring));
803  assume(l <= 0 || pLength(p) == l);
804  int i, l1;
805  poly p1 = p;
806  ring r = bucket->bucket_ring;
807
808  if (m == NULL || p == NULL) return;
809
810  if (l <= 0)
811  {
812    l1 = pLength(p1);
813    l = l1;
814  }
815  else
816    l1 = l;
817
818  kBucketMergeLm(bucket);
819  kbTest(bucket);
820  i = pLogLength(l1);
821  #ifdef USE_COEF_BUCKETS
822  number n=n_Init(1,r->cf);
823  #endif
824  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
825  {
826    //if (FALSE){
827    #ifdef USE_COEF_BUCKETS
828    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
829    {
830      number orig_coef=p_GetCoeff(bucket->coef[i],r);
831      //we take ownership:
832      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
833      number add_coef=n_Copy(p_GetCoeff(m,r),r);
834      number gcd=n_Gcd(add_coef, orig_coef,r);
835
836      if (!(n_IsOne(gcd,r)))
837      {
838        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
839        number add_coef2=n_IntDiv(add_coef, gcd,r);
840        n_Delete(&orig_coef,r);
841        n_Delete(&add_coef,r);
842        orig_coef=orig_coef2;
843        add_coef=add_coef2;
844
845        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
846        n_Delete(&n,r);
847        n=gcd;
848      }
849
850      //assume(n_IsOne(n,r));
851      number backup=p_GetCoeff(m,r);
852
853      p_SetCoeff0(m,add_coef,r);
854      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
855
856      n_Delete(&orig_coef,r);
857      p_Delete(&bucket->coef[i],r);
858
859      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
860                           bucket->buckets_length[i], l1, r);
861      l1=bucket->buckets_length[i];
862      bucket->buckets[i]=NULL;
863      bucket->buckets_length[i] = 0;
864      i = pLogLength(l1);
865      assume(l1==pLength(p1));
866
867      p_SetCoeff(m,backup,r); //deletes add_coef
868    }
869    else
870    #endif
871    {
872      MULTIPLY_BUCKET(bucket,i);
873      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
874                           bucket->buckets_length[i], l1, r);
875      l1 = bucket->buckets_length[i];
876      bucket->buckets[i] = NULL;
877      bucket->buckets_length[i] = 0;
878      i = pLogLength(l1);
879    }
880  }
881
882  else
883  {
884    #ifdef USE_COEF_BUCKETS
885    number swap_n=p_GetCoeff(m,r);
886
887    assume(n_IsOne(n,r));
888    p_SetCoeff0(m,n,r);
889    n=swap_n;
890    //p_SetCoeff0(n, swap_n, r);
891    //p_GetCoeff0(n, swap_n,r);
892    #endif
893    p1 = r->p_Procs->pp_Mult_mm(p1, m, r);
894    #ifdef USE_COEF_BUCKETS
895    //m may not be changed
896    p_SetCoeff(m,n_Copy(n,r),r);
897    #endif
898  }
899
900
901  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
902  {
903    assume(i!=0);
904    #ifdef USE_COEF_BUCKETS
905    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
906    {
907      number orig_coef=p_GetCoeff(bucket->coef[i],r);
908      //we take ownership:
909      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
910      number add_coef=n_Copy(n,r);
911      number gcd=n_Gcd(add_coef, orig_coef,r);
912
913      if (!(n_IsOne(gcd,r)))
914      {
915        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
916        number add_coef2=n_IntDiv(add_coef, gcd,r);
917        n_Delete(&orig_coef,r);
918        n_Delete(&n,r);
919        n_Delete(&add_coef,r);
920        orig_coef=orig_coef2;
921        add_coef=add_coef2;
922        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
923        n=gcd;
924      }
925      //assume(n_IsOne(n,r));
926      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
927      p1=p_Mult_nn(p1,add_coef,r);
928
929      p1 = p_Add_q(p1, bucket->buckets[i],r);
930      l1=pLength(p1);
931
932      bucket->buckets[i]=NULL;
933      n_Delete(&orig_coef,r);
934      p_Delete(&bucket->coef[i],r);
935      //l1=bucket->buckets_length[i];
936      assume(l1==pLength(p1));
937    }
938    else
939    #endif
940    {
941      //don't do that, pull out gcd
942      #ifdef USE_COEF_BUCKETS
943      if(!(n_IsOne(n,r)))
944      {
945        p1=p_Mult_nn(p1, n, r);
946        n_Delete(&n,r);
947        n=n_Init(1,r);
948      }
949      #endif
950      MULTIPLY_BUCKET(bucket,i);
951      p1 = p_Add_q(p1, bucket->buckets[i],
952                 l1, bucket->buckets_length[i], r);
953      bucket->buckets[i] = NULL;
954      bucket->buckets_length[i] = 0;
955    }
956    i = pLogLength(l1);
957  }
958
959  bucket->buckets[i] = p1;
960#ifdef USE_COEF_BUCKETS
961  assume(bucket->coef[i]==NULL);
962
963  if (!(n_IsOne(n,r)))
964  {
965    bucket->coef[i]=p_NSet(n,r);
966  }
967  else
968  {
969    bucket->coef[i]=NULL;
970    n_Delete(&n,r);
971  }
972
973  if ((p1==NULL) && (bucket->coef[i]!=NULL))
974    p_Delete(&bucket->coef[i],r);
975#endif
976  bucket->buckets_length[i]=l1;
977  if (i >= bucket->buckets_used)
978    bucket->buckets_used = i;
979  else
980    kBucketAdjustBucketsUsed(bucket);
981
982  kbTest(bucket);
983}
984
985poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
986{
987  if (q == NULL) return append;
988  poly lm;
989  loop
990  {
991    lm = kBucketGetLm(bucket);
992    if (lm == NULL) return append;
993    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
994    {
995      lm = kBucketExtractLm(bucket);
996      pNext(append) = lm;
997      pIter(append);
998    }
999    else
1000    {
1001      return append;
1002    }
1003  }
1004}
1005
1006/////////////////////////////////////////////////////////////////////////////
1007//
1008// Extract all monomials from bucket with component comp
1009// Return as a polynomial *p with length *l
1010// In other words, afterwards
1011// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
1012// and components of monomials of *p are all 0
1013//
1014
1015// Hmm... for now I'm too lazy to implement those independent of currRing
1016// But better declare it extern than including polys.h
1017extern void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r);
1018
1019void kBucketTakeOutComp(kBucket_pt bucket,
1020                        long comp,
1021                        poly *r_p, int *l)
1022{
1023  poly p = NULL, q;
1024  int i, lp = 0, lq;
1025
1026#ifndef HAVE_PSEUDO_BUCKETS
1027  kBucketMergeLm(bucket);
1028  for (i=1; i<=bucket->buckets_used; i++)
1029  {
1030    if (bucket->buckets[i] != NULL)
1031    {
1032      MULTIPLY_BUCKET(bucket,i);
1033      p_TakeOutComp(&(bucket->buckets[i]), comp, &q, &lq, bucket->bucket_ring);
1034      if (q != NULL)
1035      {
1036        assume(pLength(q) == lq);
1037        bucket->buckets_length[i] -= lq;
1038        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
1039        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1040      }
1041    }
1042  }
1043  kBucketAdjustBucketsUsed(bucket);
1044#else
1045  p_TakeOutComp(&(bucket->p), comp, &p, &lp,bucket->bucket_ring);
1046  (bucket->l) -= lp;
1047#endif
1048  *r_p = p;
1049  *l = lp;
1050
1051  kbTest(bucket);
1052}
1053
1054/////////////////////////////////////////////////////////////////////////////
1055// Reduction of Bpoly with a given poly
1056//
1057
1058extern int ksCheckCoeff(number *a, number *b);
1059
1060number kBucketPolyRed(kBucket_pt bucket,
1061                      poly p1, int l1,
1062                      poly spNoether)
1063{
1064  ring r=bucket->bucket_ring;
1065  assume((!rIsPluralRing(r))||p_LmEqual(p1,kBucketGetLm(bucket), r));
1066  assume(p1 != NULL &&
1067         p_DivisibleBy(p1,  kBucketGetLm(bucket), r));
1068  assume(pLength(p1) == (int) l1);
1069
1070  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1071  BOOLEAN reset_vec=FALSE;
1072  number rn;
1073
1074  /* we shall reduce bucket=bn*lm+... by p1=an*t+a1 where t=lm(p1)
1075     and an,bn shall be defined further down only if lc(p1)!=1
1076     we already know: an|bn and t|lm */
1077  if(a1==NULL)
1078  {
1079    p_LmDelete(&lm, r);
1080    return n_Init(1,r->cf);
1081  }
1082
1083  if (! n_IsOne(pGetCoeff(p1),r->cf))
1084  {
1085    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1086//StringSetS("##### an = "); nWrite(an); PrintS(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
1087//StringSetS("##### bn = "); nWrite(bn); PrintS(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
1088    /* ksCheckCoeff: divide out gcd from an and bn: */
1089    int ct = ksCheckCoeff(&an, &bn,r->cf);
1090    /* the previous command returns ct=0 or ct=2 iff an!=1
1091       note: an is now 1 or -1 */
1092
1093    /* setup factor for p1 which cancels leading terms */
1094    p_SetCoeff(lm, bn, r);
1095    if ((ct == 0) || (ct == 2))
1096    {
1097      /* next line used to be here before but is WRONG:
1098      kBucket_Mult_n(bucket, an);
1099        its use would result in a wrong sign for the tail of bucket
1100        in the reduction */
1101
1102      /* correct factor for cancelation by changing sign if an=-1 */
1103      if (rField_is_Ring(r))
1104        lm = p_Mult_nn(lm, an, r);
1105      else
1106        kBucket_Mult_n(bucket, an);
1107    }
1108    rn = an;
1109  }
1110  else
1111  {
1112    rn = n_Init(1,r->cf);
1113  }
1114
1115  if (p_GetComp(p1, r) != p_GetComp(lm, r))
1116  {
1117    p_SetCompP(a1, p_GetComp(lm, r), r);
1118    reset_vec = TRUE;
1119    p_SetComp(lm, p_GetComp(p1, r), r);
1120    p_Setm(lm, r);
1121  }
1122
1123  p_ExpVectorSub(lm, p1, r);
1124  l1--;
1125
1126  assume(l1==pLength(a1));
1127#if 0
1128  BOOLEAN backuped=FALSE;
1129  number coef;
1130  //@Viktor, don't ignore coefficients on monomials
1131  if(l1==1) {
1132
1133    //if (rField_is_Q(r)) {
1134      //avoid this for function fields, as gcds are expensive at the moment
1135
1136
1137      coef=p_GetCoeff(a1,r);
1138      lm=p_Mult_nn(lm, coef, r);
1139      p_SetCoeff0(a1, n_Init(1,r), r);
1140      backuped=TRUE;
1141      //WARNING: not thread_safe
1142    //deletes coef as side effect
1143    //}
1144  }
1145#endif
1146
1147  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1148
1149#if 0
1150  if (backuped)
1151    p_SetCoeff0(a1,coef,r);
1152#endif
1153
1154  p_LmDelete(&lm, r);
1155  if (reset_vec) p_SetCompP(a1, 0, r);
1156  kbTest(bucket);
1157  return rn;
1158}
1159static BOOLEAN nIsPseudoUnit(number n, ring r)
1160{
1161  if (rField_is_Zp(r))
1162    return TRUE;
1163
1164  if (rParameter(r)==NULL)
1165  {
1166    return (n_Size(n,r->cf)==1);
1167  }
1168  //if (r->parameter!=NULL)
1169  return (n_IsOne(n,r->cf) || n_IsMOne(n,r->cf));
1170}
1171
1172#ifndef USE_COEF_BUCKETS
1173void kBucketSimpleContent(kBucket_pt) {}
1174#else
1175void kBucketSimpleContent(kBucket_pt bucket)
1176{
1177  ring r=bucket->bucket_ring;
1178  int i;
1179  //PrintS("HHHHHHHHHHHHH");
1180  for (i=0;i<=MAX_BUCKET;i++)
1181  {
1182    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1183    //    PrintS("H2H2H2");
1184    if (i==0)
1185    {
1186      assume(bucket->buckets[i]==NULL);
1187    }
1188    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1189      return;
1190  }
1191  for (i=0;i<=MAX_BUCKET;i++)
1192  {
1193    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1194    //    PrintS("H2H2H2");
1195    if (i==0)
1196    {
1197      assume(bucket->buckets[i]==NULL);
1198    }
1199    if ((bucket->buckets[i]!=NULL)
1200    && (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1201      return;
1202  }
1203  //return;
1204
1205  number coef=n_Init(0,r);
1206  //ATTENTION: will not work correct for GB over ring
1207  //if (TEST_OPT_PROT)
1208  //    PrintS("CCCCCCCCCCCCC");
1209  for (i=MAX_BUCKET;i>=0;i--)
1210  {
1211    if (i==0)
1212    {
1213      assume(bucket->buckets[i]==NULL);
1214    }
1215    if (bucket->buckets[i]!=NULL)
1216    {
1217      assume(bucket->coef[i]!=NULL);
1218      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1219
1220      //in this way it should crash on programming errors, yeah
1221      number temp=n_Gcd(coef, pGetCoeff(bucket->coef[i]),r);
1222      n_Delete(&coef,r );
1223      coef=temp;
1224      if (nIsPseudoUnit(coef,r))
1225      {
1226        n_Delete(&coef,r);
1227        return;
1228      }
1229      assume(!(n_IsZero(coef,r)));
1230    }
1231  }
1232  if (n_IsZero(coef,r))
1233  {
1234    n_Delete(&coef,r);
1235    return;
1236  }
1237  if (TEST_OPT_PROT)
1238    PrintS("S");
1239  for(i=0;i<=MAX_BUCKET;i++)
1240  {
1241    if (bucket->buckets[i]!=NULL)
1242    {
1243      assume(!(n_IsZero(coef,r)));
1244      assume(bucket->coef[i]!=NULL);
1245      number lc=p_GetCoeff(bucket->coef[i],r);
1246      p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1247      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1248    }
1249  }
1250  n_Delete(&coef,r);
1251}
1252#endif
1253
1254
1255poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1256{
1257  assume(bucket->buckets[i]!=NULL);
1258
1259  poly p=bucket->buckets[i];
1260  bucket->buckets_length[i]--;
1261#ifdef USE_COEF_BUCKETS
1262  ring r=bucket->bucket_ring;
1263  if (bucket->coef[i]!=NULL)
1264  {
1265    poly next=pNext(p);
1266    if (next==NULL)
1267    {
1268      MULTIPLY_BUCKET(bucket,i);
1269      p=bucket->buckets[i];
1270      bucket->buckets[i]=NULL;
1271      return p;
1272    }
1273    else
1274    {
1275      bucket->buckets[i]=next;
1276      number c=p_GetCoeff(bucket->coef[i],r);
1277      pNext(p)=NULL;
1278      p=p_Mult_nn(p,c,r);
1279      assume(p!=NULL);
1280      return p;
1281    }
1282  }
1283  else
1284#endif
1285  {
1286    bucket->buckets[i]=pNext(bucket->buckets[i]);
1287    pNext(p)=NULL;
1288    assume(p!=NULL);
1289    return p;
1290  }
1291}
1292
1293/*
1294* input - output: a, b
1295* returns:
1296*   a := a/gcd(a,b), b := b/gcd(a,b)
1297*   and return value
1298*       0  ->  a != 1,  b != 1
1299*       1  ->  a == 1,  b != 1
1300*       2  ->  a != 1,  b == 1
1301*       3  ->  a == 1,  b == 1
1302*   this value is used to control the spolys
1303*/
1304int ksCheckCoeff(number *a, number *b, const coeffs r)
1305{
1306  int c = 0;
1307  number an = *a, bn = *b;
1308  n_Test(an,r);
1309  n_Test(bn,r);
1310
1311  number cn = n_Gcd(an, bn, r);
1312
1313  if(n_IsOne(cn, r))
1314  {
1315    an = n_Copy(an, r);
1316    bn = n_Copy(bn, r);
1317  }
1318  else
1319  {
1320    an = n_IntDiv(an, cn, r);
1321    bn = n_IntDiv(bn, cn, r);
1322  }
1323  n_Delete(&cn, r);
1324  if (n_IsOne(an, r))
1325  {
1326    c = 1;
1327  }
1328  if (n_IsOne(bn, r))
1329  {
1330    c += 2;
1331  }
1332  *a = an;
1333  *b = bn;
1334  return c;
1335}
1336
Note: See TracBrowser for help on using the repository browser.