source: git/libpolys/polys/kbuckets.cc @ 161e20

jengelh-datetimespielwiese
Last change on this file since 161e20 was 161e20, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: for PLURAL: len(p) != len(m*p)
  • Property mode set to 100644
File size: 32.1 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 ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
719  {
720    assume(pLength(bucket->buckets[i])==bucket->buckets_length[i]);
721//#ifdef USE_COEF_BUCKETS
722//     if(bucket->coef[i]!=NULL)
723//     {
724//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
725//       bucket->coef[i]=NULL;
726//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
727//                               bucket->buckets_length[i], l1,
728//                             spNoether, r);
729//     }
730//     else
731//#endif
732    MULTIPLY_BUCKET(bucket,i);
733    p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
734                            bucket->buckets_length[i], l1,
735                            spNoether, r);
736    l1 = bucket->buckets_length[i];
737    bucket->buckets[i] = NULL;
738    bucket->buckets_length[i] = 0;
739#ifdef HAVE_RINGS
740    if (rField_is_Ring(r) && !(rField_is_Domain(r)))
741    {
742      l1 = pLength(p1);
743      assume(pLength(p1) == l1);
744    }
745#endif
746#ifdef HAVE_PLURAL
747    if (rIsPluralRing(r))
748    {
749      l1 = pLength(p1);
750      assume(pLength(p1) == l1);
751    }
752#endif
753    i = pLogLength(l1);
754  }
755  else
756  {
757    pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
758    if (spNoether != NULL)
759    {
760      l1 = -1;
761      p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r);
762      i = pLogLength(l1);
763    }
764    else
765    {
766      p1 = r->p_Procs->pp_Mult_mm(p1, m, r);
767#ifdef HAVE_RINGS
768      if (rField_is_Ring(r) && !(rField_is_Domain(r)))
769      {
770        l1 = pLength(p1);
771        i = pLogLength(l1);
772      }
773#endif
774#ifdef HAVE_PLURAL
775      if (rIsPluralRing(r))
776      {
777        l1 = pLength(p1);
778        i = pLogLength(l1);
779      }
780#endif
781    }
782    pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
783  }
784
785  while (bucket->buckets[i] != NULL)
786  {
787    //kbTest(bucket);
788    MULTIPLY_BUCKET(bucket,i);
789    p1 = p_Add_q(p1, bucket->buckets[i],
790                 l1, bucket->buckets_length[i], r);
791    bucket->buckets[i] = NULL;
792    bucket->buckets_length[i] = 0;
793    i = pLogLength(l1);
794  }
795
796  bucket->buckets[i] = p1;
797  bucket->buckets_length[i]=l1;
798  if (i >= bucket->buckets_used)
799    bucket->buckets_used = i;
800  else
801    kBucketAdjustBucketsUsed(bucket);
802#else // HAVE_PSEUDO_BUCKETS
803  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
804                               bucket->l, l1,
805                               spNoether, r);
806#endif
807}
808
809//////////////////////////////////////////////////////////////////////////
810///
811/// Bpoly == Bpoly + m*p; where m is a monom
812/// Does not destroy p and m
813/// assume (l <= 0 || pLength(p) == l)
814void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
815{
816    assume((!rIsPluralRing(bucket->bucket_ring))||p_IsConstant(m, bucket->bucket_ring));
817  assume(l <= 0 || pLength(p) == l);
818  int i, l1;
819  poly p1 = p;
820  ring r = bucket->bucket_ring;
821
822  if (m == NULL || p == NULL) return;
823
824  if (l <= 0)
825  {
826    l1 = pLength(p1);
827    l = l1;
828  }
829  else
830    l1 = l;
831
832  kBucketMergeLm(bucket);
833  kbTest(bucket);
834  i = pLogLength(l1);
835  #ifdef USE_COEF_BUCKETS
836  number n=n_Init(1,r->cf);
837  #endif
838  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
839  {
840    //if (FALSE){
841    #ifdef USE_COEF_BUCKETS
842    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
843    {
844      number orig_coef=p_GetCoeff(bucket->coef[i],r);
845      //we take ownership:
846      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
847      number add_coef=n_Copy(p_GetCoeff(m,r),r);
848      number gcd=n_Gcd(add_coef, orig_coef,r);
849
850      if (!(n_IsOne(gcd,r)))
851      {
852        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
853        number add_coef2=n_IntDiv(add_coef, gcd,r);
854        n_Delete(&orig_coef,r);
855        n_Delete(&add_coef,r);
856        orig_coef=orig_coef2;
857        add_coef=add_coef2;
858
859        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
860        n_Delete(&n,r);
861        n=gcd;
862      }
863
864      //assume(n_IsOne(n,r));
865      number backup=p_GetCoeff(m,r);
866
867      p_SetCoeff0(m,add_coef,r);
868      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
869
870      n_Delete(&orig_coef,r);
871      p_Delete(&bucket->coef[i],r);
872
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      assume(l1==pLength(p1));
880
881      p_SetCoeff(m,backup,r); //deletes add_coef
882    }
883    else
884    #endif
885    {
886      MULTIPLY_BUCKET(bucket,i);
887      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
888                           bucket->buckets_length[i], l1, r);
889      l1 = bucket->buckets_length[i];
890      bucket->buckets[i] = NULL;
891      bucket->buckets_length[i] = 0;
892      i = pLogLength(l1);
893    }
894  }
895
896  else
897  {
898    #ifdef USE_COEF_BUCKETS
899    number swap_n=p_GetCoeff(m,r);
900
901    assume(n_IsOne(n,r));
902    p_SetCoeff0(m,n,r);
903    n=swap_n;
904    //p_SetCoeff0(n, swap_n, r);
905    //p_GetCoeff0(n, swap_n,r);
906    #endif
907    p1 = r->p_Procs->pp_Mult_mm(p1, m, r);
908    #ifdef USE_COEF_BUCKETS
909    //m may not be changed
910    p_SetCoeff(m,n_Copy(n,r),r);
911    #endif
912  }
913
914
915  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
916  {
917    assume(i!=0);
918    #ifdef USE_COEF_BUCKETS
919    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
920    {
921      number orig_coef=p_GetCoeff(bucket->coef[i],r);
922      //we take ownership:
923      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
924      number add_coef=n_Copy(n,r);
925      number gcd=n_Gcd(add_coef, orig_coef,r);
926
927      if (!(n_IsOne(gcd,r)))
928      {
929        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
930        number add_coef2=n_IntDiv(add_coef, gcd,r);
931        n_Delete(&orig_coef,r);
932        n_Delete(&n,r);
933        n_Delete(&add_coef,r);
934        orig_coef=orig_coef2;
935        add_coef=add_coef2;
936        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
937        n=gcd;
938      }
939      //assume(n_IsOne(n,r));
940      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
941      p1=p_Mult_nn(p1,add_coef,r);
942
943      p1 = p_Add_q(p1, bucket->buckets[i],r);
944      l1=pLength(p1);
945
946      bucket->buckets[i]=NULL;
947      n_Delete(&orig_coef,r);
948      p_Delete(&bucket->coef[i],r);
949      //l1=bucket->buckets_length[i];
950      assume(l1==pLength(p1));
951    }
952    else
953    #endif
954    {
955      //don't do that, pull out gcd
956      #ifdef USE_COEF_BUCKETS
957      if(!(n_IsOne(n,r)))
958      {
959        p1=p_Mult_nn(p1, n, r);
960        n_Delete(&n,r);
961        n=n_Init(1,r);
962      }
963      #endif
964      MULTIPLY_BUCKET(bucket,i);
965      p1 = p_Add_q(p1, bucket->buckets[i],
966                 l1, bucket->buckets_length[i], r);
967      bucket->buckets[i] = NULL;
968      bucket->buckets_length[i] = 0;
969    }
970    i = pLogLength(l1);
971  }
972
973  bucket->buckets[i] = p1;
974#ifdef USE_COEF_BUCKETS
975  assume(bucket->coef[i]==NULL);
976
977  if (!(n_IsOne(n,r)))
978  {
979    bucket->coef[i]=p_NSet(n,r);
980  }
981  else
982  {
983    bucket->coef[i]=NULL;
984    n_Delete(&n,r);
985  }
986
987  if ((p1==NULL) && (bucket->coef[i]!=NULL))
988    p_Delete(&bucket->coef[i],r);
989#endif
990  bucket->buckets_length[i]=l1;
991  if (i >= bucket->buckets_used)
992    bucket->buckets_used = i;
993  else
994    kBucketAdjustBucketsUsed(bucket);
995
996  kbTest(bucket);
997}
998
999poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
1000{
1001  if (q == NULL) return append;
1002  poly lm;
1003  loop
1004  {
1005    lm = kBucketGetLm(bucket);
1006    if (lm == NULL) return append;
1007    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
1008    {
1009      lm = kBucketExtractLm(bucket);
1010      pNext(append) = lm;
1011      pIter(append);
1012    }
1013    else
1014    {
1015      return append;
1016    }
1017  }
1018}
1019
1020/////////////////////////////////////////////////////////////////////////////
1021//
1022// Extract all monomials from bucket with component comp
1023// Return as a polynomial *p with length *l
1024// In other words, afterwards
1025// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
1026// and components of monomials of *p are all 0
1027//
1028
1029// Hmm... for now I'm too lazy to implement those independent of currRing
1030// But better declare it extern than including polys.h
1031extern void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r);
1032
1033void kBucketTakeOutComp(kBucket_pt bucket,
1034                        long comp,
1035                        poly *r_p, int *l)
1036{
1037  poly p = NULL, q;
1038  int i, lp = 0, lq;
1039
1040#ifndef HAVE_PSEUDO_BUCKETS
1041  kBucketMergeLm(bucket);
1042  for (i=1; i<=bucket->buckets_used; i++)
1043  {
1044    if (bucket->buckets[i] != NULL)
1045    {
1046      MULTIPLY_BUCKET(bucket,i);
1047      p_TakeOutComp(&(bucket->buckets[i]), comp, &q, &lq, bucket->bucket_ring);
1048      if (q != NULL)
1049      {
1050        assume(pLength(q) == lq);
1051        bucket->buckets_length[i] -= lq;
1052        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
1053        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1054      }
1055    }
1056  }
1057  kBucketAdjustBucketsUsed(bucket);
1058#else
1059  p_TakeOutComp(&(bucket->p), comp, &p, &lp,bucket->bucket_ring);
1060  (bucket->l) -= lp;
1061#endif
1062  *r_p = p;
1063  *l = lp;
1064
1065  kbTest(bucket);
1066}
1067
1068/////////////////////////////////////////////////////////////////////////////
1069// Reduction of Bpoly with a given poly
1070//
1071
1072extern int ksCheckCoeff(number *a, number *b);
1073
1074number kBucketPolyRed(kBucket_pt bucket,
1075                      poly p1, int l1,
1076                      poly spNoether)
1077{
1078  ring r=bucket->bucket_ring;
1079  assume((!rIsPluralRing(r))||p_LmEqual(p1,kBucketGetLm(bucket), r));
1080  assume(p1 != NULL &&
1081         p_DivisibleBy(p1,  kBucketGetLm(bucket), r));
1082  assume(pLength(p1) == (int) l1);
1083
1084  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1085  BOOLEAN reset_vec=FALSE;
1086  number rn;
1087
1088  /* we shall reduce bucket=bn*lm+... by p1=an*t+a1 where t=lm(p1)
1089     and an,bn shall be defined further down only if lc(p1)!=1
1090     we already know: an|bn and t|lm */
1091  if(a1==NULL)
1092  {
1093    p_LmDelete(&lm, r);
1094    return n_Init(1,r->cf);
1095  }
1096
1097  if (! n_IsOne(pGetCoeff(p1),r->cf))
1098  {
1099    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1100//StringSetS("##### an = "); nWrite(an); PrintS(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
1101//StringSetS("##### bn = "); nWrite(bn); PrintS(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
1102    /* ksCheckCoeff: divide out gcd from an and bn: */
1103    int ct = ksCheckCoeff(&an, &bn,r->cf);
1104    /* the previous command returns ct=0 or ct=2 iff an!=1
1105       note: an is now 1 or -1 */
1106
1107    /* setup factor for p1 which cancels leading terms */
1108    p_SetCoeff(lm, bn, r);
1109    if ((ct == 0) || (ct == 2))
1110    {
1111      /* next line used to be here before but is WRONG:
1112      kBucket_Mult_n(bucket, an);
1113        its use would result in a wrong sign for the tail of bucket
1114        in the reduction */
1115
1116      /* correct factor for cancelation by changing sign if an=-1 */
1117      if (rField_is_Ring(r))
1118        lm = p_Mult_nn(lm, an, r);
1119      else
1120        kBucket_Mult_n(bucket, an);
1121    }
1122    rn = an;
1123  }
1124  else
1125  {
1126    rn = n_Init(1,r->cf);
1127  }
1128
1129  if (p_GetComp(p1, r) != p_GetComp(lm, r))
1130  {
1131    p_SetCompP(a1, p_GetComp(lm, r), r);
1132    reset_vec = TRUE;
1133    p_SetComp(lm, p_GetComp(p1, r), r);
1134    p_Setm(lm, r);
1135  }
1136
1137  p_ExpVectorSub(lm, p1, r);
1138  l1--;
1139
1140  assume(l1==pLength(a1));
1141#if 0
1142  BOOLEAN backuped=FALSE;
1143  number coef;
1144  //@Viktor, don't ignore coefficients on monomials
1145  if(l1==1) {
1146
1147    //if (rField_is_Q(r)) {
1148      //avoid this for function fields, as gcds are expensive at the moment
1149
1150
1151      coef=p_GetCoeff(a1,r);
1152      lm=p_Mult_nn(lm, coef, r);
1153      p_SetCoeff0(a1, n_Init(1,r), r);
1154      backuped=TRUE;
1155      //WARNING: not thread_safe
1156    //deletes coef as side effect
1157    //}
1158  }
1159#endif
1160
1161  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1162
1163#if 0
1164  if (backuped)
1165    p_SetCoeff0(a1,coef,r);
1166#endif
1167
1168  p_LmDelete(&lm, r);
1169  if (reset_vec) p_SetCompP(a1, 0, r);
1170  kbTest(bucket);
1171  return rn;
1172}
1173static BOOLEAN nIsPseudoUnit(number n, ring r)
1174{
1175  if (rField_is_Zp(r))
1176    return TRUE;
1177
1178  if (rParameter(r)==NULL)
1179  {
1180    return (n_Size(n,r->cf)==1);
1181  }
1182  //if (r->parameter!=NULL)
1183  return (n_IsOne(n,r->cf) || n_IsMOne(n,r->cf));
1184}
1185
1186#ifndef USE_COEF_BUCKETS
1187void kBucketSimpleContent(kBucket_pt) {}
1188#else
1189void kBucketSimpleContent(kBucket_pt bucket)
1190{
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=n_Gcd(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  {
1248    n_Delete(&coef,r);
1249    return;
1250  }
1251  if (TEST_OPT_PROT)
1252    PrintS("S");
1253  for(i=0;i<=MAX_BUCKET;i++)
1254  {
1255    if (bucket->buckets[i]!=NULL)
1256    {
1257      assume(!(n_IsZero(coef,r)));
1258      assume(bucket->coef[i]!=NULL);
1259      number lc=p_GetCoeff(bucket->coef[i],r);
1260      p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1261      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1262    }
1263  }
1264  n_Delete(&coef,r);
1265}
1266#endif
1267
1268
1269poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1270{
1271  assume(bucket->buckets[i]!=NULL);
1272
1273  poly p=bucket->buckets[i];
1274  bucket->buckets_length[i]--;
1275#ifdef USE_COEF_BUCKETS
1276  ring r=bucket->bucket_ring;
1277  if (bucket->coef[i]!=NULL)
1278  {
1279    poly next=pNext(p);
1280    if (next==NULL)
1281    {
1282      MULTIPLY_BUCKET(bucket,i);
1283      p=bucket->buckets[i];
1284      bucket->buckets[i]=NULL;
1285      return p;
1286    }
1287    else
1288    {
1289      bucket->buckets[i]=next;
1290      number c=p_GetCoeff(bucket->coef[i],r);
1291      pNext(p)=NULL;
1292      p=p_Mult_nn(p,c,r);
1293      assume(p!=NULL);
1294      return p;
1295    }
1296  }
1297  else
1298#endif
1299  {
1300    bucket->buckets[i]=pNext(bucket->buckets[i]);
1301    pNext(p)=NULL;
1302    assume(p!=NULL);
1303    return p;
1304  }
1305}
1306
1307/*
1308* input - output: a, b
1309* returns:
1310*   a := a/gcd(a,b), b := b/gcd(a,b)
1311*   and return value
1312*       0  ->  a != 1,  b != 1
1313*       1  ->  a == 1,  b != 1
1314*       2  ->  a != 1,  b == 1
1315*       3  ->  a == 1,  b == 1
1316*   this value is used to control the spolys
1317*/
1318int ksCheckCoeff(number *a, number *b, const coeffs r)
1319{
1320  int c = 0;
1321  number an = *a, bn = *b;
1322  n_Test(an,r);
1323  n_Test(bn,r);
1324
1325  number cn = n_Gcd(an, bn, r);
1326
1327  if(n_IsOne(cn, r))
1328  {
1329    an = n_Copy(an, r);
1330    bn = n_Copy(bn, r);
1331  }
1332  else
1333  {
1334    an = n_IntDiv(an, cn, r);
1335    bn = n_IntDiv(bn, cn, r);
1336  }
1337  n_Delete(&cn, r);
1338  if (n_IsOne(an, r))
1339  {
1340    c = 1;
1341  }
1342  if (n_IsOne(bn, r))
1343  {
1344    c += 2;
1345  }
1346  *a = an;
1347  *b = bn;
1348  return c;
1349}
1350
Note: See TracBrowser for help on using the repository browser.