source: git/libpolys/polys/kbuckets.cc @ 803cb2d

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