source: git/libpolys/polys/kbuckets.cc @ 8e2c77

spielwiese
Last change on this file since 8e2c77 was 8e2c77, checked in by Hans Schoenemann <hannes@…>, 18 months ago
fix: redMoraNF/redMoraNFRing: n_DivBy nad SubringGCD
  • Property mode set to 100644
File size: 36.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5#include "misc/auxiliary.h"
6#include "misc/options.h"
7
8#include "polys/monomials/p_polys.h"
9#include "coeffs/coeffs.h"
10#include "coeffs/numbers.h"
11#include "polys/monomials/ring.h"
12#include "polys/kbuckets.h"
13
14#ifdef HAVE_SHIFTBBA
15#include "polys/shiftop.h"
16#endif
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],B->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(B->bucket_ring)) B->buckets_length[i] = pLength(B->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_VAR omBin kBucket_bin = omGetSpecBin(sizeof(kBucket));
46#ifdef USE_COEF_BUCKETS
47STATIC_VAR int coef_start=1;
48#endif
49//////////////////////////////////////////////////////////////////////////
50///
51/// Some internal stuff
52///
53
54// https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog
55#ifndef BUCKET_TWO_BASE
56static inline int LOG4(int v)
57{
58  const unsigned int b[] = {0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000};
59  const unsigned int S[] = {1, 2, 4, 8, 16};
60
61  unsigned int r = 0; // result of log4(v) will go here
62  if (v & b[4]) { v >>= S[4]; r |= S[3]; }
63  if (v & b[3]) { v >>= S[3]; r |= S[2]; }
64  if (v & b[2]) { v >>= S[2]; r |= S[1]; }
65  if (v & b[1]) { v >>= S[1]; r |= S[0]; }
66  return (int)r;
67}
68#endif
69
70// returns ceil(log_4(l))
71static inline unsigned int pLogLength(unsigned int l)
72{
73  unsigned int i = 0;
74
75  if (l == 0) return 0;
76  l--;
77#ifdef BUCKET_TWO_BASE
78  i=SI_LOG2(l);
79#else
80  i=LOG4(l);
81#endif
82  return i+1;
83}
84
85// returns ceil(log_4(pLength(p)))
86static inline unsigned int pLogLength(poly p)
87{
88  return pLogLength((unsigned int) pLength(p));
89}
90
91#ifdef KDEBUG
92
93#ifndef HAVE_PSEUDO_BUCKETS
94BOOLEAN kbTest_i(kBucket_pt bucket, int i)
95{//sBucketSortMerge
96  #ifdef USE_COEF_BUCKETS
97  assume(bucket->coef[0]==NULL);
98  if ((bucket->coef[i]!=NULL) && (bucket->buckets[i]==NULL))
99  {
100    dReportError("Bucket %d coef not NULL", i);
101  }
102  if (bucket->coef[i]!=NULL)
103  {
104    assume(bucket->buckets[i]!=NULL);
105    p_Test(bucket->coef[i],bucket->bucket_ring);
106  }
107  #endif
108  pFalseReturn(p_Test(bucket->buckets[i], bucket->bucket_ring));
109  if ((unsigned)bucket->buckets_length[i] != pLength(bucket->buckets[i]))
110  {
111    dReportError("Bucket %d lengths difference should:%d has:%d",
112                 i, bucket->buckets_length[i], pLength(bucket->buckets[i]));
113  }
114  else if (i > 0 && (int) pLogLength(bucket->buckets_length[i]) > i)
115  {
116    dReportError("Bucket %d too long %d",
117                 i, bucket->buckets_length[i]);
118  }
119  if (i==0 && bucket->buckets_length[0] > 1)
120  {
121    dReportError("Bucket 0 too long");
122  }
123  return TRUE;
124}
125
126
127BOOLEAN kbTest(kBucket_pt bucket)
128{
129  #ifdef HAVE_COEF_BUCKETS
130  assume(bucket->coef[0]==NULL);
131  #endif
132  int i;
133  poly lm = bucket->buckets[0];
134
135  omCheckAddrBin(bucket, kBucket_bin);
136  assume(bucket->buckets_used <= MAX_BUCKET);
137  if (! kbTest_i(bucket, 0)) return FALSE;
138  for (i=1; i<= (int) bucket->buckets_used; i++)
139  {
140    if (!kbTest_i(bucket, i)) return FALSE;
141    if (lm != NULL &&  bucket->buckets[i] != NULL
142        && p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) != 1)
143    {
144      dReportError("Bucket %d larger or equal than lm", i);
145      if (p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) ==0)
146        dReportError("Bucket %d equal to lm", i);
147      return FALSE;
148    }
149    if (!p_Test(bucket->buckets[i],bucket->bucket_ring))
150    {
151      dReportError("Bucket %d is not =0(4)", i);
152      return FALSE;
153    }
154  }
155
156  for (; i<=MAX_BUCKET; i++)
157  {
158    if (bucket->buckets[i] != NULL || bucket->buckets_length[i] != 0)
159    {
160      dReportError("Bucket %d not zero", i);
161      return FALSE;
162    }
163  }
164  for(i=0;i<=MAX_BUCKET;i++)
165  {
166    if (bucket->buckets[i]!=NULL)
167    {
168      int j;
169      for(j=i+1;j<=MAX_BUCKET;j++)
170      {
171        if (bucket->buckets[j]==bucket->buckets[i])
172        {
173          dReportError("Bucket %d %d equal", i,j);
174          return FALSE;
175        }
176      }
177    }
178    #ifdef HAVE_COEF_BUCKETS
179    if (bucket->coef[i]!=NULL)
180    {
181      int j;
182      for(j=i+1;j<=MAX_BUCKET;j++)
183      {
184        if (bucket->coef[j]==bucket->coef[i])
185        {
186          dReportError("internal coef %d %d equal", i,j);
187          return FALSE;
188        }
189      }
190    }
191    #endif
192  }
193  return TRUE;
194}
195
196#else // HAVE_PSEUDO_BUCKETS
197BOOLEAN kbTest(kBucket_pt bucket)
198{
199  return TRUE;
200}
201#endif // ! HAVE_PSEUDO_BUCKETS
202#endif // KDEBUG
203
204//////////////////////////////////////////////////////////////////////////
205///
206/// Creation/Destruction of buckets
207///
208
209kBucket_pt kBucketCreate(const ring bucket_ring)
210{
211  assume(bucket_ring != NULL);
212  kBucket_pt bucket = (kBucket_pt) omAlloc0Bin(kBucket_bin);
213  bucket->bucket_ring = bucket_ring;
214  return bucket;
215}
216void kBucketDestroy(kBucket_pt *bucket_pt)
217{
218  omFreeBin(*bucket_pt, kBucket_bin);
219  *bucket_pt = NULL;
220}
221
222
223void kBucketDeleteAndDestroy(kBucket_pt *bucket_pt)
224{
225  kBucket_pt bucket = *bucket_pt;
226  kbTest(bucket);
227  int i;
228  for (i=0; i<= bucket->buckets_used; i++)
229  {
230    p_Delete(&(bucket->buckets[i]), bucket->bucket_ring);
231#ifdef USE_COEF_BUCKETS
232    p_Delete(&(bucket->coef[i]), bucket->bucket_ring);
233#endif
234  }
235  omFreeBin(bucket, kBucket_bin);
236  *bucket_pt = NULL;
237}
238
239/////////////////////////////////////////////////////////////////////////////
240// Convertion from/to Bpolys
241//
242#ifndef HAVE_PSEUDO_BUCKETS
243
244inline void kBucketMergeLm(kBucket_pt bucket)
245{
246  kbTest(bucket);
247  if (bucket->buckets[0] != NULL)
248  {
249    poly lm = bucket->buckets[0];
250    int i = 1;
251#ifdef BUCKET_TWO_BASE
252    int l = 2;
253    while ( bucket->buckets_length[i] >= l)
254    {
255      i++;
256      l = l << 1;
257    }
258#else
259    int l = 4;
260    while ( bucket->buckets_length[i] >= l)
261    {
262      i++;
263      l = l << 2;
264    }
265#endif
266#ifndef USE_COEF_BUCKETS
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    if (i > bucket->buckets_used)  bucket->buckets_used = i;
273    bucket->buckets[0] = NULL;
274    bucket->buckets_length[0] = 0;
275    kbTest(bucket);
276#else
277    if (i > bucket->buckets_used)  bucket->buckets_used = i;
278    assume(i!=0);
279    if (bucket->buckets[i]!=NULL)
280    {
281       MULTIPLY_BUCKET(bucket,i);
282       pNext(lm) = bucket->buckets[i];
283       bucket->buckets[i] = lm;
284       bucket->buckets_length[i]++;
285       assume(i <= bucket->buckets_used+1);
286    }
287    else
288    {
289      #if 1
290       assume(bucket->buckets[i]==NULL);
291       assume(bucket->coef[0]==NULL);
292       assume(pLength(lm)==1);
293       assume(pNext(lm)==NULL);
294       number coef=p_GetCoeff(lm,bucket->bucket_ring);
295       //WARNING: not thread_safe
296       p_SetCoeff0(lm, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
297       bucket->buckets[i]=lm;
298       bucket->buckets_length[i]=1;
299       bucket->coef[i]=p_NSet(n_Copy(coef,bucket->bucket_ring),bucket->bucket_ring);
300
301       bucket->buckets[i]=lm;
302       bucket->buckets_length[i]=1;
303      #else
304       MULTIPLY_BUCKET(bucket,i);
305       pNext(lm) = bucket->buckets[i];
306       bucket->buckets[i] = lm;
307       bucket->buckets_length[i]++;
308       assume(i <= bucket->buckets_used+1);
309      #endif
310    }
311    bucket->buckets[0]=NULL;
312    bucket->buckets_length[0] = 0;
313    bucket->coef[0]=NULL;
314    kbTest(bucket);
315    #endif
316  }
317
318}
319
320BOOLEAN kBucketIsCleared(kBucket_pt bucket)
321{
322  int i;
323
324  for (i = 0;i<=MAX_BUCKET;i++)
325  {
326    if (bucket->buckets[i] != NULL) return FALSE;
327    #ifdef HAVE_COEF_BUCKETS
328    if (bucket->coef[i] != NULL) return FALSE;
329    #endif
330    if (bucket->buckets_length[i] != 0) return FALSE;
331  }
332  return TRUE;
333}
334
335void kBucketInit(kBucket_pt bucket, poly lm, int length)
336{
337  //assume(false);
338  assume(bucket != NULL);
339  assume(length <= 0 || (unsigned)length == pLength(lm));
340  assume(kBucketIsCleared(bucket));
341
342  if (lm == NULL) return;
343
344  if (length <= 0)
345    length = pLength(lm);
346
347  bucket->buckets[0] = lm;
348  #ifdef HAVE_COEF_BUCKETS
349  assume(bucket->coef[0]==NULL);
350  #endif
351  #ifdef USE_COEF_BUCKETS
352  bucket->coef[0]=NULL;
353  #endif
354  if (lm!=NULL)
355    bucket->buckets_length[0] = 1;
356  else
357    bucket->buckets_length[0]= 0;
358  if (length > 1)
359  {
360    unsigned int i = pLogLength(length-1);
361    bucket->buckets[i] = pNext(lm);
362    pNext(lm) = NULL;
363    bucket->buckets_length[i] = length-1;
364    bucket->buckets_used = i;
365  }
366  else
367  {
368    bucket->buckets_used = 0;
369  }
370}
371
372int kBucketCanonicalize(kBucket_pt bucket)
373{
374#ifndef HAVE_PSEUDO_BUCKETS
375  assume(bucket->buckets_used<=MAX_BUCKET);
376  MULTIPLY_BUCKET(bucket,1);
377  kbTest(bucket);
378  poly p = bucket->buckets[1];
379  poly lm;
380  int pl = bucket->buckets_length[1];//, i;
381  int i;
382  bucket->buckets[1] = NULL;
383  bucket->buckets_length[1] = 0;
384  #ifdef USE_COEF_BUCKETS
385    assume(bucket->coef[1]==NULL);
386  #endif
387  ring r=bucket->bucket_ring;
388
389
390  for (i=1; i<=bucket->buckets_used; i++)
391  {
392  #ifdef USE_COEF_BUCKETS
393    if (bucket->coef[i]!=NULL)
394    {
395      assume(bucket->buckets[i]!=NULL);
396      p = p_Plus_mm_Mult_qq(p, bucket->coef[i], bucket->buckets[i],
397                 pl, bucket->buckets_length[i], r);
398      p_Delete(&bucket->coef[i],r);
399      p_Delete(&bucket->buckets[i],r);
400    }
401    else
402    p = p_Add_q(p, bucket->buckets[i],
403                 pl, bucket->buckets_length[i], r);
404  #else
405    p = p_Add_q(p, bucket->buckets[i],
406                 pl, bucket->buckets_length[i], r);
407  #endif
408    if (i==1) continue;
409    bucket->buckets[i] = NULL;
410    bucket->buckets_length[i] = 0;
411  }
412  #ifdef HAVE_COEF_BUCKETS
413  assume(bucket->coef[0]==NULL);
414  #endif
415  lm = bucket->buckets[0];
416  if (lm != NULL)
417  {
418    pNext(lm) = p;
419    p = lm;
420    pl++;
421    bucket->buckets[0] = NULL;
422    bucket->buckets_length[0] = 0;
423  }
424  if (pl > 0)
425  {
426    i = pLogLength(pl);
427    bucket->buckets[i] = p;
428    bucket->buckets_length[i] = pl;
429  }
430  else
431  {
432    i = 0;
433  }
434  bucket->buckets_used = i;
435  assume(bucket->buckets_used <= MAX_BUCKET);
436  #ifdef USE_COEF_BUCKETS
437    assume(bucket->coef[0]==NULL);
438    assume(bucket->coef[i]==NULL);
439  #endif
440  assume(pLength(p) == (unsigned)pl);
441  //if (TEST_OPT_PROT) { Print("C(%d)",pl); }
442  kbTest(bucket);
443  return i;
444#endif
445}
446
447void kBucketNormalize(kBucket_pt bucket)
448{
449#ifdef HAVE_PSEUDO_BUCKETS
450  p_Normalize(bucket->p,bucket->bucket_ring);
451#else
452  MULTIPLY_BUCKET(bucket,1);
453  for (int i=0; i<=bucket->buckets_used; i++)
454  {
455    p_Normalize(bucket->buckets[i],bucket->bucket_ring);
456  }
457#endif
458}
459
460void kBucketClear(kBucket_pt bucket, poly *p, int *length)
461{
462  int i = kBucketCanonicalize(bucket);
463  if (i > 0)
464  {
465  #ifdef USE_COEF_BUCKETS
466    MULTIPLY_BUCKET(bucket,i);
467    //bucket->coef[i]=NULL;
468  #endif
469    *p = bucket->buckets[i];
470    *length = bucket->buckets_length[i];
471    bucket->buckets[i] = NULL;
472    bucket->buckets_length[i] = 0;
473    bucket->buckets_used = 0;
474
475  }
476  else
477  {
478    *p = NULL;
479    *length = 0;
480  }
481}
482
483void kBucketSetLm(kBucket_pt bucket, poly lm)
484{
485  kBucketMergeLm(bucket);
486  pNext(lm) = NULL;
487  bucket->buckets[0] = lm;
488  bucket->buckets_length[0] = 1;
489}
490
491#else // HAVE_PSEUDO_BUCKETS
492
493void kBucketInit(kBucket_pt bucket, poly lm, int length)
494{
495  int i;
496
497  assume(bucket != NULL);
498  assume(length <= 0 || length == pLength(lm));
499
500  bucket->p = lm;
501  if (length <= 0) bucket->l = pLength(lm);
502  else bucket->l = length;
503
504}
505
506const poly kBucketGetLm(kBucket_pt bucket)
507{
508  return bucket->p;
509}
510
511poly kBucketExtractLm(kBucket_pt bucket)
512{
513  poly lm = bucket->p;
514  assume(pLength(bucket->p) == bucket->l);
515  pIter(bucket->p);
516  (bucket->l)--;
517  pNext(lm) = NULL;
518  return lm;
519}
520
521void kBucketClear(kBucket_pt bucket, poly *p, int *length)
522{
523  assume(pLength(bucket->p) == bucket->l);
524  *p = bucket->p;
525  *length = bucket->l;
526  bucket->p = NULL;
527  bucket->l = 0;
528}
529
530#endif // ! HAVE_PSEUDO_BUCKETS
531//////////////////////////////////////////////////////////////////////////
532///
533/// For changing the ring of the Bpoly to new_tailBin
534///
535void kBucketShallowCopyDelete(kBucket_pt bucket,
536                              ring new_tailRing, omBin new_tailBin,
537                              pShallowCopyDeleteProc p_shallow_copy_delete)
538{
539#ifndef HAVE_PSEUDO_BUCKETS
540  int i;
541
542  kBucketCanonicalize(bucket);
543  for (i=0; i<= bucket->buckets_used; i++)
544    if (bucket->buckets[i] != NULL)
545    {
546      MULTIPLY_BUCKET(bucket,i);
547      bucket->buckets[i] = p_shallow_copy_delete(bucket->buckets[i],
548                                                 bucket->bucket_ring,
549                                                 new_tailRing,
550                                                 new_tailBin);
551    }
552#else
553  bucket->p = p_shallow_copy_delete(p,
554                                    bucket_ring,
555                                    new_tailRing,
556                                    new_tailBin);
557#endif
558  bucket->bucket_ring = new_tailRing;
559}
560
561//////////////////////////////////////////////////////////////////////////
562///
563/// Bucket number i from bucket is out of length sync, resync
564///
565void kBucketAdjust(kBucket_pt bucket, int i) {
566
567  MULTIPLY_BUCKET(bucket,i);
568
569  int l1 = bucket->buckets_length[i];
570  poly p1 = bucket->buckets[i];
571  bucket->buckets[i] = NULL;
572  bucket->buckets_length[i] = 0;
573  i = pLogLength(l1);
574
575  while (bucket->buckets[i] != NULL)
576  {
577    //kbTest(bucket);
578    MULTIPLY_BUCKET(bucket,i);
579    p1 = p_Add_q(p1, bucket->buckets[i],
580                 l1, bucket->buckets_length[i], bucket->bucket_ring);
581    bucket->buckets[i] = NULL;
582    bucket->buckets_length[i] = 0;
583    i = pLogLength(l1);
584  }
585
586  bucket->buckets[i] = p1;
587  bucket->buckets_length[i]=l1;
588  if (i >= bucket->buckets_used)
589    bucket->buckets_used = i;
590  else
591    kBucketAdjustBucketsUsed(bucket);
592}
593
594//////////////////////////////////////////////////////////////////////////
595///
596/// Multiply Bucket by number ,i.e. Bpoly == n*Bpoly
597///
598void kBucket_Mult_n(kBucket_pt bucket, number n)
599{
600#ifndef HAVE_PSEUDO_BUCKETS
601  kbTest(bucket);
602  ring r=bucket->bucket_ring;
603  int i;
604
605  for (i=0; i<= bucket->buckets_used; i++)
606  {
607    if (bucket->buckets[i] != NULL)
608    {
609#ifdef USE_COEF_BUCKETS
610      if (i<coef_start)
611        bucket->buckets[i] = __p_Mult_nn(bucket->buckets[i], n, r);
612        /* Frank Seelisch on March 11, 2010:
613           This looks a bit strange: The following "if" is indented
614           like the previous line of code. But coded as it is,
615           it should actually be two spaces less indented.
616           Question: Should the following "if" also only be
617           performed when "(i<coef_start)" is true?
618           For the time being, I leave it as it is. */
619        if (rField_is_Ring(r) && !(rField_is_Domain(r)))
620        {
621          bucket->buckets_length[i] = pLength(bucket->buckets[i]);
622          kBucketAdjust(bucket, i);
623        }
624      else
625      if (bucket->coef[i]!=NULL)
626      {
627        bucket->coef[i] = __p_Mult_nn(bucket->coef[i],n,r);
628      }
629      else
630      {
631        bucket->coef[i] = p_NSet(n_Copy(n,r),r);
632      }
633#else
634      bucket->buckets[i] = __p_Mult_nn(bucket->buckets[i], n, r);
635#endif
636    }
637  }
638  if (rField_is_Ring(r) && !(rField_is_Domain(r)))
639  {
640    for (i=0; i<= bucket->buckets_used; i++)
641    {
642      if (bucket->buckets[i] != NULL)
643      {
644        bucket->buckets_length[i] = pLength(bucket->buckets[i]);
645        kBucketAdjust(bucket, i);
646      }
647    }
648  }
649  kbTest(bucket);
650#else
651  bucket->p = __p_Mult_nn(bucket->p, n, bucket->bucket_ring);
652#endif
653}
654
655
656//////////////////////////////////////////////////////////////////////////
657///
658/// Add to Bucket a poly ,i.e. Bpoly == q+Bpoly
659///
660void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
661{
662  if (q == NULL) return;
663  assume(*l <= 0 || pLength(q) == *l);
664
665  int i, l1;
666  ring r = bucket->bucket_ring;
667
668  if (*l <= 0)
669  {
670    l1 = pLength(q);
671    *l = l1;
672  }
673  else
674    l1 = *l;
675
676  kBucketMergeLm(bucket);
677  kbTest(bucket);
678  i = pLogLength(l1);
679
680  while (bucket->buckets[i] != NULL)
681  {
682    //MULTIPLY_BUCKET(bucket,i);
683  #ifdef USE_COEF_BUCKETS
684    if (bucket->coef[i]!=NULL)
685    {
686      q = p_Plus_mm_Mult_qq(q, bucket->coef[i], bucket->buckets[i],
687                 l1, bucket->buckets_length[i], r);
688      p_Delete(&bucket->coef[i],r);
689      p_Delete(&bucket->buckets[i],r);
690    }
691    else
692    q = p_Add_q(q, bucket->buckets[i],
693                 l1, bucket->buckets_length[i], r);
694  #else
695    q = p_Add_q(q, bucket->buckets[i],
696                 l1, bucket->buckets_length[i], r);
697  #endif
698    bucket->buckets[i] = NULL;
699    bucket->buckets_length[i] = 0;
700    i = pLogLength(l1);
701    assume(i<= MAX_BUCKET);
702    assume(bucket->buckets_used<= MAX_BUCKET);
703  }
704
705  kbTest(bucket);
706  bucket->buckets[i] = q;
707  bucket->buckets_length[i]=l1;
708  if (i >= bucket->buckets_used)
709    bucket->buckets_used = i;
710  else
711    kBucketAdjustBucketsUsed(bucket);
712  kbTest(bucket);
713}
714
715
716
717//////////////////////////////////////////////////////////////////////////
718///
719/// Bpoly == Bpoly - m*p; where m is a monom
720/// Does not destroy p and m
721/// assume (*l <= 0 || pLength(p) == *l)
722void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l,
723                            poly spNoether)
724{
725  assume(*l <= 0 || pLength(p) == *l);
726  int i, l1;
727  poly p1 = p;
728  ring r = bucket->bucket_ring;
729
730  if (*l <= 0)
731  {
732    l1 = pLength(p1);
733    *l = l1;
734  }
735  else
736    l1 = *l;
737
738  if (m == NULL || p == NULL) return;
739
740#ifndef HAVE_PSEUDO_BUCKETS
741  kBucketMergeLm(bucket);
742  kbTest(bucket);
743  i = pLogLength(l1);
744
745#if defined(HAVE_PLURAL)
746  if ((rField_is_Ring(r) && !(rField_is_Domain(r)))
747  ||(rIsPluralRing(r)))
748  {
749    pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
750    p1=r->p_Procs->pp_mm_Mult(p,m,r);
751    pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
752    l1=pLength(p1);
753    i = pLogLength(l1);
754  }
755  else
756#endif
757  {
758    if ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
759    {
760      assume(pLength(bucket->buckets[i])==(unsigned)bucket->buckets_length[i]);
761//#ifdef USE_COEF_BUCKETS
762//     if(bucket->coef[i]!=NULL)
763//     {
764//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
765//       bucket->coef[i]=NULL;
766//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
767//                               bucket->buckets_length[i], l1,
768//                             spNoether, r);
769//     }
770//     else
771//#endif
772      MULTIPLY_BUCKET(bucket,i);
773      p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
774                            bucket->buckets_length[i], l1,
775                            spNoether, r);
776      l1 = bucket->buckets_length[i];
777      bucket->buckets[i] = NULL;
778      bucket->buckets_length[i] = 0;
779      i = pLogLength(l1);
780    }
781    else
782    {
783      pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
784      if (spNoether != NULL)
785      {
786        l1 = -1;
787        p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r);
788        i = pLogLength(l1);
789      }
790      else
791      {
792        p1 = r->p_Procs->pp_mm_Mult(p1, m, r);
793      }
794      pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
795    }
796  }
797
798  while (bucket->buckets[i] != NULL)
799  {
800    //kbTest(bucket);
801    MULTIPLY_BUCKET(bucket,i);
802    p1 = p_Add_q(p1, bucket->buckets[i],
803                 l1, bucket->buckets_length[i], r);
804    bucket->buckets[i] = NULL;
805    bucket->buckets_length[i] = 0;
806    i = pLogLength(l1);
807  }
808
809  bucket->buckets[i] = p1;
810  bucket->buckets_length[i]=l1;
811  if (i >= bucket->buckets_used)
812    bucket->buckets_used = i;
813  else
814    kBucketAdjustBucketsUsed(bucket);
815#else // HAVE_PSEUDO_BUCKETS
816  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
817                               bucket->l, l1,
818                               spNoether, r);
819#endif
820}
821
822//////////////////////////////////////////////////////////////////////////
823///
824/// Bpoly == Bpoly + m*p; where m is a monom
825/// Does not destroy p and m
826/// assume (l <= 0 || pLength(p) == l)
827void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
828{
829    assume((!rIsPluralRing(bucket->bucket_ring))||p_IsConstant(m, bucket->bucket_ring));
830  assume(l <= 0 || pLength(p) == (unsigned)l);
831  int i, l1;
832  poly p1 = p;
833  ring r = bucket->bucket_ring;
834
835  if (m == NULL || p == NULL) return;
836
837  if (l <= 0)
838  {
839    l1 = pLength(p1);
840    l = l1;
841  }
842  else
843    l1 = l;
844
845  kBucketMergeLm(bucket);
846  kbTest(bucket);
847  i = pLogLength(l1);
848  #ifdef USE_COEF_BUCKETS
849  number n=n_Init(1,r->cf);
850  #endif
851  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
852  {
853    //if (FALSE){
854    #ifdef USE_COEF_BUCKETS
855    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
856    {
857      number orig_coef=p_GetCoeff(bucket->coef[i],r);
858      //we take ownership:
859      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
860      number add_coef=n_Copy(p_GetCoeff(m,r),r);
861      number gcd=n_Gcd(add_coef, orig_coef,r);
862
863      if (!(n_IsOne(gcd,r)))
864      {
865        number orig_coef2=n_ExactDiv(orig_coef,gcd,r);
866        number add_coef2=n_ExactDiv(add_coef, gcd,r);
867        n_Delete(&orig_coef,r);
868        n_Delete(&add_coef,r);
869        orig_coef=orig_coef2;
870        add_coef=add_coef2;
871
872        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
873        n_Delete(&n,r);
874        n=gcd;
875      }
876
877      //assume(n_IsOne(n,r));
878      number backup=p_GetCoeff(m,r);
879
880      p_SetCoeff0(m,add_coef,r);
881      bucket->buckets[i]=__p_Mult_nn(bucket->buckets[i],orig_coef,r);
882
883      n_Delete(&orig_coef,r);
884      p_Delete(&bucket->coef[i],r);
885
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      assume(l1==pLength(p1));
893
894      p_SetCoeff(m,backup,r); //deletes add_coef
895    }
896    else
897    #endif
898    {
899      MULTIPLY_BUCKET(bucket,i);
900      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
901                           bucket->buckets_length[i], l1, r);
902      l1 = bucket->buckets_length[i];
903      bucket->buckets[i] = NULL;
904      bucket->buckets_length[i] = 0;
905      i = pLogLength(l1);
906    }
907  }
908  else
909  {
910    #ifdef USE_COEF_BUCKETS
911    number swap_n=p_GetCoeff(m,r);
912
913    assume(n_IsOne(n,r));
914    p_SetCoeff0(m,n,r);
915    n=swap_n;
916    //p_SetCoeff0(n, swap_n, r);
917    //p_GetCoeff0(n, swap_n,r);
918    #endif
919    p1 = r->p_Procs->pp_Mult_mm(p1, m, r);
920    #ifdef USE_COEF_BUCKETS
921    //m may not be changed
922    p_SetCoeff(m,n_Copy(n,r),r);
923    #endif
924  }
925
926  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
927  {
928    assume(i!=0);
929    #ifdef USE_COEF_BUCKETS
930    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
931    {
932      number orig_coef=p_GetCoeff(bucket->coef[i],r);
933      //we take ownership:
934      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
935      number add_coef=n_Copy(n,r);
936      number gcd=n_Gcd(add_coef, orig_coef,r);
937
938      if (!(n_IsOne(gcd,r)))
939      {
940        number orig_coef2=n_ExactDiv(orig_coef,gcd,r);
941        number add_coef2=n_ExactDiv(add_coef, gcd,r);
942        n_Delete(&orig_coef,r);
943        n_Delete(&n,r);
944        n_Delete(&add_coef,r);
945        orig_coef=orig_coef2;
946        add_coef=add_coef2;
947        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
948        n=gcd;
949      }
950      //assume(n_IsOne(n,r));
951      bucket->buckets[i]=__p_Mult_nn(bucket->buckets[i],orig_coef,r);
952      p1=__p_Mult_nn(p1,add_coef,r);
953
954      p1 = p_Add_q(p1, bucket->buckets[i],r);
955      l1=pLength(p1);
956
957      bucket->buckets[i]=NULL;
958      n_Delete(&orig_coef,r);
959      p_Delete(&bucket->coef[i],r);
960      //l1=bucket->buckets_length[i];
961      assume(l1==pLength(p1));
962    }
963    else
964    #endif
965    {
966      //don't do that, pull out gcd
967      #ifdef USE_COEF_BUCKETS
968      if(!(n_IsOne(n,r)))
969      {
970        p1=__p_Mult_nn(p1, n, r);
971        n_Delete(&n,r);
972        n=n_Init(1,r);
973      }
974      #endif
975      MULTIPLY_BUCKET(bucket,i);
976      p1 = p_Add_q(p1, bucket->buckets[i],
977                 l1, bucket->buckets_length[i], r);
978      bucket->buckets[i] = NULL;
979      bucket->buckets_length[i] = 0;
980    }
981    i = pLogLength(l1);
982  }
983
984  bucket->buckets[i] = p1;
985#ifdef USE_COEF_BUCKETS
986  assume(bucket->coef[i]==NULL);
987
988  if (!(n_IsOne(n,r)))
989  {
990    bucket->coef[i]=p_NSet(n,r);
991  }
992  else
993  {
994    bucket->coef[i]=NULL;
995    n_Delete(&n,r);
996  }
997
998  if (p1==NULL)
999    p_Delete(&bucket->coef[i],r);
1000#endif
1001  bucket->buckets_length[i]=l1;
1002  if (i > bucket->buckets_used)
1003    bucket->buckets_used = i;
1004  else
1005    kBucketAdjustBucketsUsed(bucket);
1006
1007  kbTest(bucket);
1008}
1009
1010poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
1011{
1012  if (q == NULL) return append;
1013  poly lm;
1014  loop
1015  {
1016    lm = kBucketGetLm(bucket);
1017    if (lm == NULL) return append;
1018    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
1019    {
1020      lm = kBucketExtractLm(bucket);
1021      pNext(append) = lm;
1022      pIter(append);
1023    }
1024    else
1025    {
1026      return append;
1027    }
1028  }
1029}
1030
1031/////////////////////////////////////////////////////////////////////////////
1032//
1033// Extract all monomials from bucket with component comp
1034// Return as a polynomial *p with length *l
1035// In other words, afterwards
1036// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
1037// and components of monomials of *p are all 0
1038//
1039
1040// Hmm... for now I'm too lazy to implement those independent of currRing
1041// But better declare it extern than including polys.h
1042extern void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r);
1043
1044void kBucketTakeOutComp(kBucket_pt bucket,
1045                        long comp,
1046                        poly *r_p, int *l)
1047{
1048  poly p = NULL, q;
1049  int i, lp = 0, lq;
1050
1051#ifndef HAVE_PSEUDO_BUCKETS
1052  kBucketMergeLm(bucket);
1053  for (i=1; i<=bucket->buckets_used; i++)
1054  {
1055    if (bucket->buckets[i] != NULL)
1056    {
1057      MULTIPLY_BUCKET(bucket,i);
1058      p_TakeOutComp(&(bucket->buckets[i]), comp, &q, &lq, bucket->bucket_ring);
1059      if (q != NULL)
1060      {
1061        assume(pLength(q) == (unsigned)lq);
1062        bucket->buckets_length[i] -= lq;
1063        assume(pLength(bucket->buckets[i]) == (unsigned)bucket->buckets_length[i]);
1064        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1065      }
1066    }
1067  }
1068  kBucketAdjustBucketsUsed(bucket);
1069#else
1070  p_TakeOutComp(&(bucket->p), comp, &p, &lp,bucket->bucket_ring);
1071  (bucket->l) -= lp;
1072#endif
1073  *r_p = p;
1074  *l = lp;
1075
1076  kbTest(bucket);
1077}
1078
1079/////////////////////////////////////////////////////////////////////////////
1080// Reduction of Bpoly with a given poly
1081//
1082
1083number kBucketPolyRed(kBucket_pt bucket,
1084                      poly p1, int l1,
1085                      poly spNoether)
1086{
1087  ring r=bucket->bucket_ring;
1088  assume((!rIsPluralRing(r))||p_LmEqual(p1,kBucketGetLm(bucket), r));
1089  assume(p1 != NULL &&
1090         p_DivisibleBy(p1,  kBucketGetLm(bucket), r));
1091  assume(pLength(p1) == (unsigned) l1);
1092
1093  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1094  BOOLEAN reset_vec=FALSE;
1095  number rn;
1096
1097  /* we shall reduce bucket=bn*lm+... by p1=an*t+a1 where t=lm(p1)
1098     and an,bn shall be defined further down only if lc(p1)!=1
1099     we already know: an|bn and t|lm */
1100  if(a1==NULL)
1101  {
1102    p_LmDelete(&lm, r);
1103    return n_Init(1,r->cf);
1104  }
1105
1106  if (! n_IsOne(pGetCoeff(p1),r->cf))
1107  {
1108    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1109    /* ksCheckCoeff: divide out gcd from an and bn: */
1110    int ct = ksCheckCoeff(&an, &bn,r->cf);
1111    /* the previous command returns ct=0 or ct=2 iff an!=1
1112       note: an is now 1 or -1 */
1113
1114    /* setup factor for p1 which cancels leading terms */
1115    p_SetCoeff(lm, bn, r);
1116    if ((ct == 0) || (ct == 2))
1117    {
1118      /* correct factor for cancelation by changing sign if an=-1 */
1119      if (rField_is_Ring(r))
1120        lm = __p_Mult_nn(lm, an, r);
1121      else
1122        kBucket_Mult_n(bucket, an);
1123    }
1124    rn = an;
1125  }
1126  else
1127  {
1128    rn = n_Init(1,r->cf);
1129  }
1130
1131  if (p_GetComp(p1, r) != p_GetComp(lm, r))
1132  {
1133    p_SetCompP(a1, p_GetComp(lm, r), r);
1134    reset_vec = TRUE;
1135    p_SetComp(lm, p_GetComp(p1, r), r);
1136    p_Setm(lm, r);
1137  }
1138
1139  p_ExpVectorSub(lm, p1, r);
1140  l1--;
1141
1142  assume((unsigned)l1==pLength(a1));
1143
1144#ifdef HAVE_SHIFTBBA
1145  poly lmRight;
1146  poly lm_org;
1147  if (r->isLPring)
1148  {
1149    int firstBlock = p_mFirstVblock(p1, r);
1150    lm_org=lm;
1151    k_SplitFrame(lm, lmRight, si_max(firstBlock, 1), r);
1152  }
1153#endif
1154#if 0
1155  BOOLEAN backuped=FALSE;
1156  number coef;
1157  //@Viktor, don't ignore coefficients on monomials
1158  if(l1==1) {
1159
1160    //if (rField_is_Q(r)) {
1161      //avoid this for function fields, as gcds are expensive at the moment
1162
1163
1164      coef=p_GetCoeff(a1,r);
1165      lm=p_Mult_nn(lm, coef, r);
1166      p_SetCoeff0(a1, n_Init(1,r), r);
1167      backuped=TRUE;
1168      //WARNING: not thread_safe
1169    //deletes coef as side effect
1170    //}
1171  }
1172#endif
1173
1174#ifdef HAVE_SHIFTBBA
1175  if (r->isLPring)
1176  {
1177    poly tmp=r->p_Procs->pp_Mult_mm(a1, lmRight, r);
1178    kBucket_Minus_m_Mult_p(bucket, lm,tmp, &l1, spNoether);
1179    p_Delete(&tmp,r);
1180    p_LmDelete(&lmRight,r);
1181    p_LmDelete(lm_org,r);
1182  }
1183  else
1184#endif
1185  {
1186    kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1187  }
1188
1189#if 0
1190  if (backuped)
1191    p_SetCoeff0(a1,coef,r);
1192#endif
1193
1194  p_LmDelete(&lm, r);
1195  if (reset_vec) p_SetCompP(a1, 0, r);
1196  kbTest(bucket);
1197  return rn;
1198}
1199
1200void kBucketPolyRedNF(kBucket_pt bucket,
1201                      poly p1, int l1,
1202                      poly spNoether)
1203{
1204  ring r=bucket->bucket_ring;
1205  assume((!rIsPluralRing(r))||p_LmEqual(p1,kBucketGetLm(bucket), r));
1206  assume(p1 != NULL &&
1207         p_DivisibleBy(p1,  kBucketGetLm(bucket), r));
1208  assume(pLength(p1) == (unsigned) l1);
1209
1210  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1211  BOOLEAN reset_vec=FALSE;
1212
1213  /* we shall reduce bucket=bn*lm+... by p1=an*t+a1 where t=lm(p1)
1214     and an,bn shall be defined further down only if lc(p1)!=1
1215     we already know: an|bn and t|lm */
1216  if(a1==NULL)
1217  {
1218    p_LmDelete(&lm, r);
1219    return;
1220  }
1221
1222  #ifdef KDEBUG
1223  if (n_DivBy(pGetCoeff(lm),pGetCoeff(p1),r->cf))
1224  #endif
1225  {
1226    number c=n_Div(pGetCoeff(lm),pGetCoeff(p1),r->cf);
1227    //StringSetS("mult cf:");n_Write(c,r->cf);StringAppendS("\n");
1228    //PrintS(StringEndS());
1229    #ifdef KDEBUG
1230    if (n_IsZero(c,r->cf))
1231    {
1232      StringSetS("a/b: ");n_Write(pGetCoeff(lm),r->cf);StringAppendS(" / ");
1233      n_Write(pGetCoeff(p1),r->cf);StringAppendS("\n");PrintS(StringEndS());
1234    }
1235    #endif
1236    p_SetCoeff(lm,c,r);
1237  }
1238  #ifdef KDEBUG
1239  else
1240  {
1241    PrintS("bug\n");
1242  }
1243  #endif
1244  if (p_GetComp(p1, r) != p_GetComp(lm, r))
1245  {
1246    p_SetCompP(a1, p_GetComp(lm, r), r);
1247    reset_vec = TRUE;
1248    p_SetComp(lm, p_GetComp(p1, r), r);
1249    p_Setm(lm, r);
1250  }
1251
1252  p_ExpVectorSub(lm, p1, r);
1253  l1--;
1254
1255  assume((unsigned)l1==pLength(a1));
1256
1257#ifdef HAVE_SHIFTBBA
1258  poly lmRight;
1259  poly lm_org;
1260  if (r->isLPring)
1261  {
1262    int firstBlock = p_mFirstVblock(p1, r);
1263    lm_org=lm;
1264    k_SplitFrame(lm, lmRight, si_max(firstBlock, 1), r);
1265  }
1266#endif
1267
1268#ifdef HAVE_SHIFTBBA
1269  if (r->isLPring)
1270  {
1271    poly tmp=r->p_Procs->pp_Mult_mm(a1, lmRight, r);
1272    kBucket_Minus_m_Mult_p(bucket, lm,tmp, &l1, spNoether);
1273    p_Delete(&tmp,r);
1274    p_LmDelete(&lmRight,r);
1275    p_LmDelete(lm_org,r);
1276  }
1277  else
1278#endif
1279  {
1280    kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1281  }
1282
1283#if 0
1284  if (backuped)
1285    p_SetCoeff0(a1,coef,r);
1286#endif
1287
1288  p_LmDelete(&lm, r);
1289  if (reset_vec) p_SetCompP(a1, 0, r);
1290  kbTest(bucket);
1291  return;
1292}
1293
1294#ifndef USE_COEF_BUCKETS
1295void kBucketSimpleContent(kBucket_pt bucket)
1296{
1297  if (bucket->buckets[0]==NULL) return;
1298
1299  ring r=bucket->bucket_ring;
1300  if (rField_is_Ring(r)) return;
1301
1302  coeffs cf=r->cf;
1303  if (cf->cfSubringGcd==ndGcd) /* trivial gcd*/ return;
1304
1305  number nn=pGetCoeff(bucket->buckets[0]);
1306  //if ((bucket->buckets_used==0)
1307  //&&(!n_IsOne(nn,cf)))
1308  //{
1309  //  if (TEST_OPT_PROT) PrintS("@");
1310  //  p_SetCoeff(bucket->buckets[0],n_Init(1,cf),r);
1311  //  return;
1312  //}
1313
1314  if (n_Size(nn,cf)<2) return;
1315
1316  //kBucketAdjustBucketsUsed(bucket);
1317  number coef=n_Copy(nn,cf);
1318  // find an initial guess of a gcd
1319  for (int i=1; i<=bucket->buckets_used;i++)
1320  {
1321    if (bucket->buckets[i]!=NULL)
1322    {
1323      number t=p_InitContent(bucket->buckets[i],r);
1324      if (n_Size(t,cf)<2)
1325      {
1326        n_Delete(&t,cf);
1327        n_Delete(&coef,cf);
1328        return;
1329      }
1330      number t2=n_SubringGcd(coef,t,cf);
1331      n_Delete(&t,cf);
1332      n_Delete(&coef,cf);
1333      coef=t2;
1334      if (n_Size(coef,cf)<2) { n_Delete(&coef,cf);return;}
1335    }
1336  }
1337  // find the gcd
1338  for (int i=0; i<=bucket->buckets_used;i++)
1339  {
1340    if (bucket->buckets[i]!=NULL)
1341    {
1342      poly p=bucket->buckets[i];
1343      while(p!=NULL)
1344      {
1345        number t=n_SubringGcd(coef,pGetCoeff(p),cf);
1346        if (n_Size(t,cf)<2)
1347        {
1348          n_Delete(&t,cf);
1349          n_Delete(&coef,cf);
1350          return;
1351        }
1352        pIter(p);
1353      }
1354    }
1355  }
1356  // divided by the gcd
1357  if (TEST_OPT_PROT) PrintS("@");
1358  for (int i=bucket->buckets_used;i>=0;i--)
1359  {
1360    if (bucket->buckets[i]!=NULL)
1361    {
1362      poly p=bucket->buckets[i];
1363      while(p!=NULL)
1364      {
1365        number d = n_ExactDiv(pGetCoeff(p),coef,cf);
1366        p_SetCoeff(p,d,r);
1367        pIter(p);
1368      }
1369    }
1370  }
1371  n_Delete(&coef,cf);
1372}
1373#else
1374static BOOLEAN nIsPseudoUnit(number n, ring r)
1375{
1376  if (rField_is_Zp(r))
1377    return TRUE;
1378
1379  if (rParameter(r)==NULL)
1380  {
1381    return (n_Size(n,r->cf)==1);
1382  }
1383  //if (r->parameter!=NULL)
1384  return (n_IsOne(n,r->cf) || n_IsMOne(n,r->cf));
1385}
1386
1387void kBucketSimpleContent(kBucket_pt bucket)
1388{
1389  ring r=bucket->bucket_ring;
1390  int i;
1391  //PrintS("HHHHHHHHHHHHH");
1392  for (i=0;i<=MAX_BUCKET;i++)
1393  {
1394    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1395    //    PrintS("H2H2H2");
1396    if (i==0)
1397    {
1398      assume(bucket->buckets[i]==NULL);
1399    }
1400    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1401      return;
1402  }
1403  for (i=0;i<=MAX_BUCKET;i++)
1404  {
1405    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1406    //    PrintS("H2H2H2");
1407    if (i==0)
1408    {
1409      assume(bucket->buckets[i]==NULL);
1410    }
1411    if ((bucket->buckets[i]!=NULL)
1412    && (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1413      return;
1414  }
1415  //return;
1416
1417  number coef=n_Init(0,r);
1418  //ATTENTION: will not work correct for GB over ring
1419  //if (TEST_OPT_PROT)
1420  //    PrintS("CCCCCCCCCCCCC");
1421  for (i=MAX_BUCKET;i>=0;i--)
1422  {
1423    if (i==0)
1424    {
1425      assume(bucket->buckets[i]==NULL);
1426    }
1427    if (bucket->buckets[i]!=NULL)
1428    {
1429      assume(bucket->coef[i]!=NULL);
1430      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1431
1432      //in this way it should crash on programming errors, yeah
1433      number temp=n_Gcd(coef, pGetCoeff(bucket->coef[i]),r);
1434      n_Delete(&coef,r );
1435      coef=temp;
1436      if (nIsPseudoUnit(coef,r))
1437      {
1438        n_Delete(&coef,r);
1439        return;
1440      }
1441      assume(!(n_IsZero(coef,r)));
1442    }
1443  }
1444  if (n_IsZero(coef,r))
1445  {
1446    n_Delete(&coef,r);
1447    return;
1448  }
1449  if (TEST_OPT_PROT)
1450    PrintS("S");
1451  for(i=0;i<=MAX_BUCKET;i++)
1452  {
1453    if (bucket->buckets[i]!=NULL)
1454    {
1455      assume(!(n_IsZero(coef,r)));
1456      assume(bucket->coef[i]!=NULL);
1457      number lc=p_GetCoeff(bucket->coef[i],r);
1458      p_SetCoeff(bucket->coef[i], n_ExactDiv(lc,coef,r),r);
1459      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1460    }
1461  }
1462  n_Delete(&coef,r);
1463}
1464#endif
1465
1466
1467poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1468{
1469  assume(bucket->buckets[i]!=NULL);
1470
1471  poly p=bucket->buckets[i];
1472  bucket->buckets_length[i]--;
1473#ifdef USE_COEF_BUCKETS
1474  ring r=bucket->bucket_ring;
1475  if (bucket->coef[i]!=NULL)
1476  {
1477    poly next=pNext(p);
1478    if (next==NULL)
1479    {
1480      MULTIPLY_BUCKET(bucket,i);
1481      p=bucket->buckets[i];
1482      bucket->buckets[i]=NULL;
1483      return p;
1484    }
1485    else
1486    {
1487      bucket->buckets[i]=next;
1488      number c=p_GetCoeff(bucket->coef[i],r);
1489      pNext(p)=NULL;
1490      p=__p_Mult_nn(p,c,r);
1491      assume(p!=NULL);
1492      return p;
1493    }
1494  }
1495  else
1496#endif
1497  {
1498    bucket->buckets[i]=pNext(bucket->buckets[i]);
1499    pNext(p)=NULL;
1500    assume(p!=NULL);
1501    return p;
1502  }
1503}
1504
1505/*
1506* input - output: a, b
1507* returns:
1508*   a := a/gcd(a,b), b := b/gcd(a,b)
1509*   and return value
1510*       0  ->  a != 1,  b != 1
1511*       1  ->  a == 1,  b != 1
1512*       2  ->  a != 1,  b == 1
1513*       3  ->  a == 1,  b == 1
1514*   this value is used to control the spolys
1515*/
1516int ksCheckCoeff(number *a, number *b, const coeffs r)
1517{
1518  int c = 0;
1519  number an = *a, bn = *b;
1520  n_Test(an,r);
1521  n_Test(bn,r);
1522
1523  if (UNLIKELY(nCoeff_is_Ring(r) && n_DivBy(bn,an,r))) // in NF
1524  {
1525      bn = n_ExactDiv(bn, an, r);
1526      an = n_Init(1, r);
1527  }
1528  else
1529  {
1530    number cn = n_SubringGcd(an, bn, r);
1531    if(n_IsOne(cn, r))
1532    {
1533      an = n_Copy(an, r); // a/1
1534      bn = n_Copy(bn, r); // b/1
1535    }
1536    else
1537    {
1538      an = n_ExactDiv(an, cn, r);
1539      bn = n_ExactDiv(bn, cn, r);
1540    }
1541    n_Delete(&cn, r);
1542  }
1543  if (n_IsOne(an, r))
1544  {
1545    c = 1;
1546  }
1547  if (n_IsOne(bn, r))
1548  {
1549    c += 2;
1550  }
1551  *a = an;
1552  *b = bn;
1553  return c;
1554}
1555
Note: See TracBrowser for help on using the repository browser.