source: git/kernel/kbuckets.cc @ d772c3

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