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

spielwiese
Last change on this file since e6f1e6 was e6f1e6, checked in by Hans Schoenemann <hannes@…>, 8 years ago
code cleanup: HAVE_RINGS is mostly not needed
  • Property mode set to 100644
File size: 31.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4//#include <kernel/mod2.h>
5
6
7
8
9#include <omalloc/omalloc.h>
10#include <misc/auxiliary.h>
11
12#include <polys/monomials/p_polys.h>
13//#include <kernel/pShallowCopyDelete.h>
14#include <coeffs/coeffs.h>
15#include <polys/monomials/ring.h>
16//#include <kernel/GBEngine/kutil.h>
17#include <polys/kbuckets.h>
18
19// #include <polys/operations/pShallowCopyDelete.h>
20
21
22#ifdef HAVE_COEF_BUCKETS
23#define USE_COEF_BUCKETS
24#endif
25
26#ifdef USE_COEF_BUCKETS
27#ifdef HAVE_RINGS_OLD
28#define MULTIPLY_BUCKET(B,I) do                                        \
29  { if (B->coef[I]!=NULL)                                              \
30    {                                                                  \
31      assume(p_IsConstant(B->Coef[i],B->bucket->ring));           \
32      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
33      B->coef[I]=NULL;                                                 \
34    }                                                                  \
35  } while(0)                                                           \
36    if (rField_is_Ring(B->bucket_ring)) B->buckets_length[i] = pLength(B->buckets[i]);
37#else
38#define MULTIPLY_BUCKET(B,I) do                                        \
39  { if (B->coef[I]!=NULL)                                              \
40    {                                                                  \
41      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
42      B->coef[I]=NULL;                                                 \
43    }                                                                  \
44  } while(0)
45#endif
46#else
47#define MULTIPLY_BUCKET(B,I)
48#endif
49static omBin kBucket_bin = omGetSpecBin(sizeof(kBucket));
50#ifdef USE_COEF_BUCKETS
51static int coef_start=1;
52#endif
53//////////////////////////////////////////////////////////////////////////
54///
55/// Some internal stuff
56///
57
58// returns ceil(log_4(l))
59inline unsigned int pLogLength(unsigned int l)
60{
61  unsigned int i = 0;
62
63  if (l == 0) return 0;
64  l--;
65#ifdef BUCKET_TWO_BASE
66  while ((l = (l >> 1))) i++;
67#else
68  while ((l = (l >> 2))) i++;
69#endif
70  return i+1;
71}
72
73// returns ceil(log_4(pLength(p)))
74inline unsigned int pLogLength(poly p)
75{
76  return pLogLength((unsigned int) pLength(p));
77}
78
79#ifdef KDEBUG
80
81#ifndef HAVE_PSEUDO_BUCKETS
82BOOLEAN kbTest_i(kBucket_pt bucket, int i)
83{//sBucketSortMerge
84  #ifdef USE_COEF_BUCKETS
85  assume(bucket->coef[0]==NULL);
86  if ((bucket->coef[i]!=NULL) && (bucket->buckets[i]==NULL))
87  {
88    dReportError("Bucket %d coef not NULL", i);
89  }
90  if (bucket->coef[i]!=NULL)
91  {
92    assume(bucket->buckets[i]!=NULL);
93    p_Test(bucket->coef[i],bucket->bucket_ring);
94  }
95  #endif
96  pFalseReturn(p_Test(bucket->buckets[i], bucket->bucket_ring));
97  if (bucket->buckets_length[i] != pLength(bucket->buckets[i]))
98  {
99    dReportError("Bucket %d lengths difference should:%d has:%d",
100                 i, bucket->buckets_length[i], pLength(bucket->buckets[i]));
101  }
102  else if (i > 0 && (int) pLogLength(bucket->buckets_length[i]) > i)
103  {
104    dReportError("Bucket %d too long %d",
105                 i, bucket->buckets_length[i]);
106  }
107  if (i==0 && bucket->buckets_length[0] > 1)
108  {
109    dReportError("Bucket 0 too long");
110  }
111  return TRUE;
112}
113
114
115BOOLEAN kbTest(kBucket_pt bucket)
116{
117  #ifdef HAVE_COEF_BUCKETS
118  assume(bucket->coef[0]==NULL);
119  #endif
120  int i;
121  poly lm = bucket->buckets[0];
122
123  omCheckAddrBin(bucket, kBucket_bin);
124  assume(bucket->buckets_used <= MAX_BUCKET);
125  if (! kbTest_i(bucket, 0)) return FALSE;
126  for (i=1; i<= (int) bucket->buckets_used; i++)
127  {
128    if (!kbTest_i(bucket, i)) return FALSE;
129    if (lm != NULL &&  bucket->buckets[i] != NULL
130        && p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) != 1)
131    {
132      dReportError("Bucket %d larger or equal than lm", i);
133      if (p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) ==0)
134        dReportError("Bucket %d equal to lm", i);
135      return FALSE;
136    }
137    if (!p_Test(bucket->buckets[i],bucket->bucket_ring))
138    {
139      dReportError("Bucket %d is not =0(4)", i);
140      return FALSE;
141    }
142  }
143
144  for (; i<=MAX_BUCKET; i++)
145  {
146    if (bucket->buckets[i] != NULL || bucket->buckets_length[i] != 0)
147    {
148      dReportError("Bucket %d not zero", i);
149      return FALSE;
150    }
151  }
152  for(i=0;i<=MAX_BUCKET;i++)
153  {
154    if (bucket->buckets[i]!=NULL)
155    {
156      int j;
157      for(j=i+1;j<=MAX_BUCKET;j++)
158      {
159        if (bucket->buckets[j]==bucket->buckets[i])
160        {
161          dReportError("Bucket %d %d equal", i,j);
162          return FALSE;
163        }
164      }
165    }
166    #ifdef HAVE_COEF_BUCKETS
167    if (bucket->coef[i]!=NULL)
168    {
169      int j;
170      for(j=i+1;j<=MAX_BUCKET;j++)
171      {
172        if (bucket->coef[j]==bucket->coef[i])
173        {
174          dReportError("internal coef %d %d equal", i,j);
175          return FALSE;
176        }
177      }
178    }
179    #endif
180  }
181  return TRUE;
182}
183
184#else // HAVE_PSEUDO_BUCKETS
185BOOLEAN kbTest(kBucket_pt bucket)
186{
187  return TRUE;
188}
189#endif // ! HAVE_PSEUDO_BUCKETS
190#endif // KDEBUG
191
192//////////////////////////////////////////////////////////////////////////
193///
194/// Creation/Destruction of buckets
195///
196
197kBucket_pt kBucketCreate(const ring bucket_ring)
198{
199  assume(bucket_ring != NULL);
200  kBucket_pt bucket = (kBucket_pt) omAlloc0Bin(kBucket_bin);
201  bucket->bucket_ring = bucket_ring;
202  return bucket;
203}
204void kBucketDestroy(kBucket_pt *bucket_pt)
205{
206  omFreeBin(*bucket_pt, kBucket_bin);
207  *bucket_pt = NULL;
208}
209
210
211void kBucketDeleteAndDestroy(kBucket_pt *bucket_pt)
212{
213  kBucket_pt bucket = *bucket_pt;
214  kbTest(bucket);
215  int i;
216  for (i=0; i<= bucket->buckets_used; i++)
217  {
218
219    if (bucket->buckets[i] != NULL)
220    {
221      p_Delete(&(bucket->buckets[i]), bucket->bucket_ring);
222#ifdef USE_COEF_BUCKETS
223      if (bucket->coef[i]!=NULL)
224        p_Delete(&(bucket->coef[i]), bucket->bucket_ring);
225#endif
226    }
227  }
228  omFreeBin(bucket, kBucket_bin);
229  *bucket_pt = NULL;
230}
231
232/////////////////////////////////////////////////////////////////////////////
233// Convertion from/to Bpolys
234//
235#ifndef HAVE_PSEUDO_BUCKETS
236
237inline void kBucketMergeLm(kBucket_pt bucket)
238{
239  kbTest(bucket);
240  if (bucket->buckets[0] != NULL)
241  {
242    poly lm = bucket->buckets[0];
243    int i = 1;
244#ifdef BUCKET_TWO_BASE
245    int l = 2;
246    while ( bucket->buckets_length[i] >= l)
247    {
248      i++;
249      l = l << 1;
250    }
251#else
252    int l = 4;
253    while ( bucket->buckets_length[i] >= l)
254    {
255      i++;
256      l = l << 2;
257    }
258#endif
259#ifndef USE_COEF_BUCKETS
260    MULTIPLY_BUCKET(bucket,i);
261    pNext(lm) = bucket->buckets[i];
262    bucket->buckets[i] = lm;
263    bucket->buckets_length[i]++;
264    assume(i <= bucket->buckets_used+1);
265    if (i > bucket->buckets_used)  bucket->buckets_used = i;
266    bucket->buckets[0] = NULL;
267    bucket->buckets_length[0] = 0;
268    kbTest(bucket);
269#else
270    if (i > bucket->buckets_used)  bucket->buckets_used = i;
271    assume(i!=0);
272    if (bucket->buckets[i]!=NULL)
273    {
274       MULTIPLY_BUCKET(bucket,i);
275       pNext(lm) = bucket->buckets[i];
276       bucket->buckets[i] = lm;
277       bucket->buckets_length[i]++;
278       assume(i <= bucket->buckets_used+1);
279    }
280    else
281    {
282      #if 1
283       assume(bucket->buckets[i]==NULL);
284       assume(bucket->coef[0]==NULL);
285       assume(pLength(lm)==1);
286       assume(pNext(lm)==NULL);
287       number coef=p_GetCoeff(lm,bucket->bucket_ring);
288       //WARNING: not thread_safe
289       p_SetCoeff0(lm, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
290       bucket->buckets[i]=lm;
291       bucket->buckets_length[i]=1;
292       bucket->coef[i]=p_NSet(n_Copy(coef,bucket->bucket_ring),bucket->bucket_ring);
293
294       bucket->buckets[i]=lm;
295       bucket->buckets_length[i]=1;
296      #else
297       MULTIPLY_BUCKET(bucket,i);
298       pNext(lm) = bucket->buckets[i];
299       bucket->buckets[i] = lm;
300       bucket->buckets_length[i]++;
301       assume(i <= bucket->buckets_used+1);
302      #endif
303    }
304    bucket->buckets[0]=NULL;
305    bucket->buckets_length[0] = 0;
306    bucket->coef[0]=NULL;
307    kbTest(bucket);
308    #endif
309  }
310
311}
312
313BOOLEAN kBucketIsCleared(kBucket_pt bucket)
314{
315  int i;
316
317  for (i = 0;i<=MAX_BUCKET;i++)
318  {
319    if (bucket->buckets[i] != NULL) return FALSE;
320    #ifdef HAVE_COEF_BUCKETS
321    if (bucket->coef[i] != NULL) return FALSE;
322    #endif
323    if (bucket->buckets_length[i] != 0) return FALSE;
324  }
325  return TRUE;
326}
327
328void kBucketInit(kBucket_pt bucket, poly lm, int length)
329{
330  //assume(false);
331  assume(bucket != NULL);
332  assume(length <= 0 || length == pLength(lm));
333  assume(kBucketIsCleared(bucket));
334
335  if (lm == NULL) return;
336
337  if (length <= 0)
338    length = pLength(lm);
339
340  bucket->buckets[0] = lm;
341  #ifdef HAVE_COEF_BUCKETS
342  assume(bucket->coef[0]==NULL);
343  #endif
344  #ifdef USE_COEF_BUCKETS
345  bucket->coef[0]=NULL;
346  #endif
347  if (lm!=NULL)
348    bucket->buckets_length[0] = 1;
349  else
350    bucket->buckets_length[0]= 0;
351  if (length > 1)
352  {
353    unsigned int i = pLogLength(length-1);
354    bucket->buckets[i] = pNext(lm);
355    pNext(lm) = NULL;
356    bucket->buckets_length[i] = length-1;
357    bucket->buckets_used = i;
358  }
359  else
360  {
361    bucket->buckets_used = 0;
362  }
363}
364
365int kBucketCanonicalize(kBucket_pt bucket)
366{
367  assume(bucket->buckets_used<=MAX_BUCKET);
368  MULTIPLY_BUCKET(bucket,1);
369  kbTest(bucket);
370  poly p = bucket->buckets[1];
371  poly lm;
372  int pl = bucket->buckets_length[1];//, i;
373  int i;
374  bucket->buckets[1] = NULL;
375  bucket->buckets_length[1] = 0;
376  #ifdef USE_COEF_BUCKETS
377    assume(bucket->coef[1]==NULL);
378  #endif
379  ring r=bucket->bucket_ring;
380
381
382  for (i=1; i<=bucket->buckets_used; i++)
383  {
384  #ifdef USE_COEF_BUCKETS
385    if (bucket->coef[i]!=NULL)
386    {
387      assume(bucket->buckets[i]!=NULL);
388      p = p_Plus_mm_Mult_qq(p, bucket->coef[i], bucket->buckets[i],
389                 pl, bucket->buckets_length[i], r);
390      p_Delete(&bucket->coef[i],r);
391      p_Delete(&bucket->buckets[i],r);
392    }
393    else
394    p = p_Add_q(p, bucket->buckets[i],
395                 pl, bucket->buckets_length[i], r);
396  #else
397    p = p_Add_q(p, bucket->buckets[i],
398                 pl, bucket->buckets_length[i], r);
399  #endif
400    if (i==1) continue;
401    bucket->buckets[i] = NULL;
402    bucket->buckets_length[i] = 0;
403  }
404  #ifdef HAVE_COEF_BUCKETS
405  assume(bucket->coef[0]==NULL);
406  #endif
407  lm = bucket->buckets[0];
408  if (lm != NULL)
409  {
410    pNext(lm) = p;
411    p = lm;
412    pl++;
413    bucket->buckets[0] = NULL;
414    bucket->buckets_length[0] = 0;
415  }
416  if (pl > 0)
417  {
418    i = pLogLength(pl);
419    bucket->buckets[i] = p;
420    bucket->buckets_length[i] = pl;
421  }
422  else
423  {
424    i = 0;
425  }
426  bucket->buckets_used = i;
427  assume(bucket->buckets_used <= MAX_BUCKET);
428  #ifdef USE_COEF_BUCKETS
429    assume(bucket->coef[0]==NULL);
430    assume(bucket->coef[i]==NULL);
431  #endif
432  assume(pLength(p) == (int) pl);
433  //if (TEST_OPT_PROT) { Print("C(%d)",pl); }
434  kbTest(bucket);
435  return i;
436}
437
438void kBucketClear(kBucket_pt bucket, poly *p, int *length)
439{
440  int i = kBucketCanonicalize(bucket);
441  if (i > 0)
442  {
443  #ifdef USE_COEF_BUCKETS
444    MULTIPLY_BUCKET(bucket,i);
445    //bucket->coef[i]=NULL;
446  #endif
447    *p = bucket->buckets[i];
448    *length = bucket->buckets_length[i];
449    bucket->buckets[i] = NULL;
450    bucket->buckets_length[i] = 0;
451    bucket->buckets_used = 0;
452
453  }
454  else
455  {
456    *p = NULL;
457    *length = 0;
458  }
459}
460
461void kBucketSetLm(kBucket_pt bucket, poly lm)
462{
463  kBucketMergeLm(bucket);
464  pNext(lm) = NULL;
465  bucket->buckets[0] = lm;
466  bucket->buckets_length[0] = 1;
467}
468
469#else // HAVE_PSEUDO_BUCKETS
470
471void kBucketInit(kBucket_pt bucket, poly lm, int length)
472{
473  int i;
474
475  assume(bucket != NULL);
476  assume(length <= 0 || length == pLength(lm));
477
478  bucket->p = lm;
479  if (length <= 0) bucket->l = pLength(lm);
480  else bucket->l = length;
481
482}
483
484const poly kBucketGetLm(kBucket_pt bucket)
485{
486  return bucket->p;
487}
488
489poly kBucketExtractLm(kBucket_pt bucket)
490{
491  poly lm = bucket->p;
492  assume(pLength(bucket->p) == bucket->l);
493  pIter(bucket->p);
494  (bucket->l)--;
495  pNext(lm) = NULL;
496  return lm;
497}
498
499void kBucketClear(kBucket_pt bucket, poly *p, int *length)
500{
501  assume(pLength(bucket->p) == bucket->l);
502  *p = bucket->p;
503  *length = bucket->l;
504  bucket->p = NULL;
505  bucket->l = 0;
506}
507
508#endif // ! HAVE_PSEUDO_BUCKETS
509//////////////////////////////////////////////////////////////////////////
510///
511/// For changing the ring of the Bpoly to new_tailBin
512///
513void kBucketShallowCopyDelete(kBucket_pt bucket,
514                              ring new_tailRing, omBin new_tailBin,
515                              pShallowCopyDeleteProc p_shallow_copy_delete)
516{
517#ifndef HAVE_PSEUDO_BUCKETS
518  int i;
519
520  kBucketCanonicalize(bucket);
521  for (i=0; i<= bucket->buckets_used; i++)
522    if (bucket->buckets[i] != NULL)
523    {
524      MULTIPLY_BUCKET(bucket,i);
525      bucket->buckets[i] = p_shallow_copy_delete(bucket->buckets[i],
526                                                 bucket->bucket_ring,
527                                                 new_tailRing,
528                                                 new_tailBin);
529    }
530#else
531  bucket->p = p_shallow_copy_delete(p,
532                                    bucket_ring,
533                                    new_tailRing,
534                                    new_tailBin);
535#endif
536  bucket->bucket_ring = new_tailRing;
537}
538
539//////////////////////////////////////////////////////////////////////////
540///
541/// Bucket number i from bucket is out of length sync, resync
542///
543void kBucketAdjust(kBucket_pt bucket, int i) {
544
545  MULTIPLY_BUCKET(bucket,i);
546
547  int l1 = bucket->buckets_length[i];
548  poly p1 = bucket->buckets[i];
549  bucket->buckets[i] = NULL;
550  bucket->buckets_length[i] = 0;
551  i = pLogLength(l1);
552
553  while (bucket->buckets[i] != NULL)
554  {
555    //kbTest(bucket);
556    MULTIPLY_BUCKET(bucket,i);
557    p1 = p_Add_q(p1, bucket->buckets[i],
558                 l1, bucket->buckets_length[i], bucket->bucket_ring);
559    bucket->buckets[i] = NULL;
560    bucket->buckets_length[i] = 0;
561    i = pLogLength(l1);
562  }
563
564  bucket->buckets[i] = p1;
565  bucket->buckets_length[i]=l1;
566  if (i >= bucket->buckets_used)
567    bucket->buckets_used = i;
568  else
569    kBucketAdjustBucketsUsed(bucket);
570}
571
572//////////////////////////////////////////////////////////////////////////
573///
574/// Multiply Bucket by number ,i.e. Bpoly == n*Bpoly
575///
576void kBucket_Mult_n(kBucket_pt bucket, number n)
577{
578#ifndef HAVE_PSEUDO_BUCKETS
579  kbTest(bucket);
580  ring r=bucket->bucket_ring;
581  int i;
582
583  for (i=0; i<= bucket->buckets_used; i++)
584  {
585    if (bucket->buckets[i] != NULL)
586    {
587#ifdef USE_COEF_BUCKETS
588      if (i<coef_start)
589        bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
590        /* Frank Seelisch on March 11, 2010:
591           This looks a bit strange: The following "if" is indented
592           like the previous line of code. But coded as it is,
593           it should actually be two spaces less indented.
594           Question: Should the following "if" also only be
595           performed when "(i<coef_start)" is true?
596           For the time being, I leave it as it is. */
597        if (rField_is_Ring(r) && !(rField_is_Domain(r)))
598        {
599          bucket->buckets_length[i] = pLength(bucket->buckets[i]);
600          kBucketAdjust(bucket, i);
601        }
602      else
603      if (bucket->coef[i]!=NULL)
604      {
605        bucket->coef[i] = p_Mult_nn(bucket->coef[i],n,r);
606      }
607      else
608      {
609        bucket->coef[i] = p_NSet(n_Copy(n,r),r);
610      }
611#else
612      bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
613      if (rField_is_Ring(r) && !(rField_is_Domain(r)))
614      {
615        bucket->buckets_length[i] = pLength(bucket->buckets[i]);
616        kBucketAdjust(bucket, i);
617      }
618#endif
619    }
620  }
621  kbTest(bucket);
622#else
623  bucket->p = p_Mult_nn(bucket->p, n, bucket->bucket_ring);
624#endif
625}
626
627
628//////////////////////////////////////////////////////////////////////////
629///
630/// Add to Bucket a poly ,i.e. Bpoly == q+Bpoly
631///
632void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
633{
634  if (q == NULL) return;
635  assume(*l <= 0 || pLength(q) == *l);
636
637  int i, l1;
638  ring r = bucket->bucket_ring;
639
640  if (*l <= 0)
641  {
642    l1 = pLength(q);
643    *l = l1;
644  }
645  else
646    l1 = *l;
647
648  kBucketMergeLm(bucket);
649  kbTest(bucket);
650  i = pLogLength(l1);
651
652  while (bucket->buckets[i] != NULL)
653  {
654    //MULTIPLY_BUCKET(bucket,i);
655  #ifdef USE_COEF_BUCKETS
656    if (bucket->coef[i]!=NULL)
657    {
658      q = p_Plus_mm_Mult_qq(q, bucket->coef[i], bucket->buckets[i],
659                 l1, bucket->buckets_length[i], r);
660      p_Delete(&bucket->coef[i],r);
661      p_Delete(&bucket->buckets[i],r);
662    }
663    else
664    q = p_Add_q(q, bucket->buckets[i],
665                 l1, bucket->buckets_length[i], r);
666  #else
667    q = p_Add_q(q, bucket->buckets[i],
668                 l1, bucket->buckets_length[i], r);
669  #endif
670    bucket->buckets[i] = NULL;
671    bucket->buckets_length[i] = 0;
672    i = pLogLength(l1);
673    assume(i<= MAX_BUCKET);
674    assume(bucket->buckets_used<= MAX_BUCKET);
675  }
676
677  kbTest(bucket);
678  bucket->buckets[i] = q;
679  bucket->buckets_length[i]=l1;
680  if (i >= bucket->buckets_used)
681    bucket->buckets_used = i;
682  else
683    kBucketAdjustBucketsUsed(bucket);
684  kbTest(bucket);
685}
686
687
688
689//////////////////////////////////////////////////////////////////////////
690///
691/// Bpoly == Bpoly - m*p; where m is a monom
692/// Does not destroy p and m
693/// assume (*l <= 0 || pLength(p) == *l)
694void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l,
695                            poly spNoether)
696{
697  assume(*l <= 0 || pLength(p) == *l);
698  int i, l1;
699  poly p1 = p;
700  ring r = bucket->bucket_ring;
701
702  if (*l <= 0)
703  {
704    l1 = pLength(p1);
705    *l = l1;
706  }
707  else
708    l1 = *l;
709
710  if (m == NULL || p == NULL) return;
711
712#ifndef HAVE_PSEUDO_BUCKETS
713  kBucketMergeLm(bucket);
714  kbTest(bucket);
715  i = pLogLength(l1);
716
717#if defined(HAVE_PLURAL)
718  if ((rField_is_Ring(r) && !(rField_is_Domain(r)))
719  ||(rIsPluralRing(r)))
720  {
721    pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
722    p1=pp_Mult_mm(p,m,r);
723    pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
724    l1=pLength(p1);
725    i = pLogLength(l1);
726  }
727  else
728#endif
729  {
730    if ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
731    {
732      assume(pLength(bucket->buckets[i])==bucket->buckets_length[i]);
733//#ifdef USE_COEF_BUCKETS
734//     if(bucket->coef[i]!=NULL)
735//     {
736//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
737//       bucket->coef[i]=NULL;
738//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
739//                               bucket->buckets_length[i], l1,
740//                             spNoether, r);
741//     }
742//     else
743//#endif
744      MULTIPLY_BUCKET(bucket,i);
745      p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
746                            bucket->buckets_length[i], l1,
747                            spNoether, r);
748      l1 = bucket->buckets_length[i];
749      bucket->buckets[i] = NULL;
750      bucket->buckets_length[i] = 0;
751      i = pLogLength(l1);
752    }
753    else
754    {
755      pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
756      if (spNoether != NULL)
757      {
758        l1 = -1;
759        p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r);
760        i = pLogLength(l1);
761      }
762      else
763      {
764        p1 = r->p_Procs->pp_Mult_mm(p1, m, r);
765      }
766      pSetCoeff0(m, n_InpNeg(pGetCoeff(m),r->cf));
767    }
768  }
769
770  while (bucket->buckets[i] != NULL)
771  {
772    //kbTest(bucket);
773    MULTIPLY_BUCKET(bucket,i);
774    p1 = p_Add_q(p1, bucket->buckets[i],
775                 l1, bucket->buckets_length[i], r);
776    bucket->buckets[i] = NULL;
777    bucket->buckets_length[i] = 0;
778    i = pLogLength(l1);
779  }
780
781  bucket->buckets[i] = p1;
782  bucket->buckets_length[i]=l1;
783  if (i >= bucket->buckets_used)
784    bucket->buckets_used = i;
785  else
786    kBucketAdjustBucketsUsed(bucket);
787#else // HAVE_PSEUDO_BUCKETS
788  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
789                               bucket->l, l1,
790                               spNoether, r);
791#endif
792}
793
794//////////////////////////////////////////////////////////////////////////
795///
796/// Bpoly == Bpoly + m*p; where m is a monom
797/// Does not destroy p and m
798/// assume (l <= 0 || pLength(p) == l)
799void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
800{
801    assume((!rIsPluralRing(bucket->bucket_ring))||p_IsConstant(m, bucket->bucket_ring));
802  assume(l <= 0 || pLength(p) == l);
803  int i, l1;
804  poly p1 = p;
805  ring r = bucket->bucket_ring;
806
807  if (m == NULL || p == NULL) return;
808
809  if (l <= 0)
810  {
811    l1 = pLength(p1);
812    l = l1;
813  }
814  else
815    l1 = l;
816
817  kBucketMergeLm(bucket);
818  kbTest(bucket);
819  i = pLogLength(l1);
820  #ifdef USE_COEF_BUCKETS
821  number n=n_Init(1,r->cf);
822  #endif
823  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
824  {
825    //if (FALSE){
826    #ifdef USE_COEF_BUCKETS
827    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
828    {
829      number orig_coef=p_GetCoeff(bucket->coef[i],r);
830      //we take ownership:
831      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
832      number add_coef=n_Copy(p_GetCoeff(m,r),r);
833      number gcd=n_Gcd(add_coef, orig_coef,r);
834
835      if (!(n_IsOne(gcd,r)))
836      {
837        number orig_coef2=n_ExactDiv(orig_coef,gcd,r);
838        number add_coef2=n_ExactDiv(add_coef, gcd,r);
839        n_Delete(&orig_coef,r);
840        n_Delete(&add_coef,r);
841        orig_coef=orig_coef2;
842        add_coef=add_coef2;
843
844        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
845        n_Delete(&n,r);
846        n=gcd;
847      }
848
849      //assume(n_IsOne(n,r));
850      number backup=p_GetCoeff(m,r);
851
852      p_SetCoeff0(m,add_coef,r);
853      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
854
855      n_Delete(&orig_coef,r);
856      p_Delete(&bucket->coef[i],r);
857
858      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
859                           bucket->buckets_length[i], l1, r);
860      l1=bucket->buckets_length[i];
861      bucket->buckets[i]=NULL;
862      bucket->buckets_length[i] = 0;
863      i = pLogLength(l1);
864      assume(l1==pLength(p1));
865
866      p_SetCoeff(m,backup,r); //deletes add_coef
867    }
868    else
869    #endif
870    {
871      MULTIPLY_BUCKET(bucket,i);
872      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
873                           bucket->buckets_length[i], l1, r);
874      l1 = bucket->buckets_length[i];
875      bucket->buckets[i] = NULL;
876      bucket->buckets_length[i] = 0;
877      i = pLogLength(l1);
878    }
879  }
880  else
881  {
882    #ifdef USE_COEF_BUCKETS
883    number swap_n=p_GetCoeff(m,r);
884
885    assume(n_IsOne(n,r));
886    p_SetCoeff0(m,n,r);
887    n=swap_n;
888    //p_SetCoeff0(n, swap_n, r);
889    //p_GetCoeff0(n, swap_n,r);
890    #endif
891    p1 = r->p_Procs->pp_Mult_mm(p1, m, r);
892    #ifdef USE_COEF_BUCKETS
893    //m may not be changed
894    p_SetCoeff(m,n_Copy(n,r),r);
895    #endif
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_ExactDiv(orig_coef,gcd,r);
913        number add_coef2=n_ExactDiv(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(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
1084//StringSetS("##### bn = "); nWrite(bn); PrintS(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
1085    /* ksCheckCoeff: divide out gcd from an and bn: */
1086    int ct = ksCheckCoeff(&an, &bn,r->cf);
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#if 0
1125  BOOLEAN backuped=FALSE;
1126  number coef;
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 0
1147  if (backuped)
1148    p_SetCoeff0(a1,coef,r);
1149#endif
1150
1151  p_LmDelete(&lm, r);
1152  if (reset_vec) p_SetCompP(a1, 0, r);
1153  kbTest(bucket);
1154  return rn;
1155}
1156
1157#ifndef USE_COEF_BUCKETS
1158void kBucketSimpleContent(kBucket_pt) {}
1159#else
1160static BOOLEAN nIsPseudoUnit(number n, ring r)
1161{
1162  if (rField_is_Zp(r))
1163    return TRUE;
1164
1165  if (rParameter(r)==NULL)
1166  {
1167    return (n_Size(n,r->cf)==1);
1168  }
1169  //if (r->parameter!=NULL)
1170  return (n_IsOne(n,r->cf) || n_IsMOne(n,r->cf));
1171}
1172
1173void kBucketSimpleContent(kBucket_pt bucket)
1174{
1175  ring r=bucket->bucket_ring;
1176  int i;
1177  //PrintS("HHHHHHHHHHHHH");
1178  for (i=0;i<=MAX_BUCKET;i++)
1179  {
1180    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1181    //    PrintS("H2H2H2");
1182    if (i==0)
1183    {
1184      assume(bucket->buckets[i]==NULL);
1185    }
1186    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1187      return;
1188  }
1189  for (i=0;i<=MAX_BUCKET;i++)
1190  {
1191    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1192    //    PrintS("H2H2H2");
1193    if (i==0)
1194    {
1195      assume(bucket->buckets[i]==NULL);
1196    }
1197    if ((bucket->buckets[i]!=NULL)
1198    && (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1199      return;
1200  }
1201  //return;
1202
1203  number coef=n_Init(0,r);
1204  //ATTENTION: will not work correct for GB over ring
1205  //if (TEST_OPT_PROT)
1206  //    PrintS("CCCCCCCCCCCCC");
1207  for (i=MAX_BUCKET;i>=0;i--)
1208  {
1209    if (i==0)
1210    {
1211      assume(bucket->buckets[i]==NULL);
1212    }
1213    if (bucket->buckets[i]!=NULL)
1214    {
1215      assume(bucket->coef[i]!=NULL);
1216      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1217
1218      //in this way it should crash on programming errors, yeah
1219      number temp=n_Gcd(coef, pGetCoeff(bucket->coef[i]),r);
1220      n_Delete(&coef,r );
1221      coef=temp;
1222      if (nIsPseudoUnit(coef,r))
1223      {
1224        n_Delete(&coef,r);
1225        return;
1226      }
1227      assume(!(n_IsZero(coef,r)));
1228    }
1229  }
1230  if (n_IsZero(coef,r))
1231  {
1232    n_Delete(&coef,r);
1233    return;
1234  }
1235  if (TEST_OPT_PROT)
1236    PrintS("S");
1237  for(i=0;i<=MAX_BUCKET;i++)
1238  {
1239    if (bucket->buckets[i]!=NULL)
1240    {
1241      assume(!(n_IsZero(coef,r)));
1242      assume(bucket->coef[i]!=NULL);
1243      number lc=p_GetCoeff(bucket->coef[i],r);
1244      p_SetCoeff(bucket->coef[i], n_ExactDiv(lc,coef,r),r);
1245      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1246    }
1247  }
1248  n_Delete(&coef,r);
1249}
1250#endif
1251
1252
1253poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1254{
1255  assume(bucket->buckets[i]!=NULL);
1256
1257  poly p=bucket->buckets[i];
1258  bucket->buckets_length[i]--;
1259#ifdef USE_COEF_BUCKETS
1260  ring r=bucket->bucket_ring;
1261  if (bucket->coef[i]!=NULL)
1262  {
1263    poly next=pNext(p);
1264    if (next==NULL)
1265    {
1266      MULTIPLY_BUCKET(bucket,i);
1267      p=bucket->buckets[i];
1268      bucket->buckets[i]=NULL;
1269      return p;
1270    }
1271    else
1272    {
1273      bucket->buckets[i]=next;
1274      number c=p_GetCoeff(bucket->coef[i],r);
1275      pNext(p)=NULL;
1276      p=p_Mult_nn(p,c,r);
1277      assume(p!=NULL);
1278      return p;
1279    }
1280  }
1281  else
1282#endif
1283  {
1284    bucket->buckets[i]=pNext(bucket->buckets[i]);
1285    pNext(p)=NULL;
1286    assume(p!=NULL);
1287    return p;
1288  }
1289}
1290
1291/*
1292* input - output: a, b
1293* returns:
1294*   a := a/gcd(a,b), b := b/gcd(a,b)
1295*   and return value
1296*       0  ->  a != 1,  b != 1
1297*       1  ->  a == 1,  b != 1
1298*       2  ->  a != 1,  b == 1
1299*       3  ->  a == 1,  b == 1
1300*   this value is used to control the spolys
1301*/
1302int ksCheckCoeff(number *a, number *b, const coeffs r)
1303{
1304  int c = 0;
1305  number an = *a, bn = *b;
1306  n_Test(an,r);
1307  n_Test(bn,r);
1308
1309  number cn = n_SubringGcd(an, bn, r);
1310
1311  if(n_IsOne(cn, r))
1312  {
1313    an = n_Copy(an, r);
1314    bn = n_Copy(bn, r);
1315  }
1316  else
1317  {
1318    an = n_Div(an, cn, r); n_Normalize(an,r);
1319    bn = n_Div(bn, cn, r); n_Normalize(bn,r);
1320  }
1321  n_Delete(&cn, r);
1322  if (n_IsOne(an, r))
1323  {
1324    c = 1;
1325  }
1326  if (n_IsOne(bn, r))
1327  {
1328    c += 2;
1329  }
1330  *a = an;
1331  *b = bn;
1332  return c;
1333}
1334
Note: See TracBrowser for help on using the repository browser.