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

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