source: git/libpolys/polys/kbuckets.cc @ 975db18

spielwiese
Last change on this file since 975db18 was a3d94c, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: finxing a circular dependency due to kbuckets.h and pShallowCopyDelete.h
  • Property mode set to 100644
File size: 31.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6//#include <kernel/mod2.h>
7#include <omalloc/omalloc.h>
8#include <polys/monomials/p_polys.h>
9//#include <kernel/febase.h>
10//#include <kernel/pShallowCopyDelete.h>
11#include <coeffs/coeffs.h>
12#include <polys/monomials/ring.h>
13//#include <kernel/p_Procs.h>
14//#include <kernel/kutil.h>
15#include <polys/kbuckets.h>
16
17// #include <polys/operations/pShallowCopyDelete.h>
18
19
20#ifdef HAVE_COEF_BUCKETS
21#define USE_COEF_BUCKETS
22#endif
23
24#ifdef USE_COEF_BUCKETS
25#ifdef HAVE_RINGS_OLD
26#define MULTIPLY_BUCKET(B,I) do                                        \
27  { if (B->coef[I]!=NULL)                                              \
28    {                                                                  \
29      assume(p_IsConstant(B->Coef[i],B->bucket->ring));           \
30      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
31      B->coef[I]=NULL;                                                 \
32    }                                                                  \
33  } while(0)                                                           \
34    if (rField_is_Ring(B->bucket_ring)) B->buckets_length[i] = pLength(B->buckets[i]);
35#else
36#define MULTIPLY_BUCKET(B,I) do                                        \
37  { if (B->coef[I]!=NULL)                                              \
38    {                                                                  \
39      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
40      B->coef[I]=NULL;                                                 \
41    }                                                                  \
42  } while(0)
43#endif
44#else
45#define MULTIPLY_BUCKET(B,I)
46#endif
47static omBin kBucket_bin = omGetSpecBin(sizeof(kBucket));
48#ifdef USE_COEF_BUCKETS
49static int coef_start=1;
50#endif
51//////////////////////////////////////////////////////////////////////////
52///
53/// Some internal stuff
54///
55
56// returns ceil(log_4(l))
57inline unsigned int pLogLength(unsigned int l)
58{
59  unsigned int i = 0;
60
61  if (l == 0) return 0;
62  l--;
63#ifdef BUCKET_TWO_BASE
64  while ((l = (l >> 1))) i++;
65#else
66  while ((l = (l >> 2))) i++;
67#endif
68  return i+1;
69}
70
71// returns ceil(log_4(pLength(p)))
72inline unsigned int pLogLength(poly p)
73{
74  return pLogLength((unsigned int) pLength(p));
75}
76
77#ifdef KDEBUG
78
79#ifndef HAVE_PSEUDO_BUCKETS
80BOOLEAN kbTest_i(kBucket_pt bucket, int i)
81{//sBucketSortMerge
82  #ifdef USE_COEF_BUCKETS
83  assume(bucket->coef[0]==NULL);
84  if ((bucket->coef[i]!=NULL) && (bucket->buckets[i]==NULL))
85  {
86    dReportError("Bucket %d coef not NULL", i);
87  }
88  if (bucket->coef[i]!=NULL)
89  {
90    assume(bucket->buckets[i]!=NULL);
91    _p_Test(bucket->coef[i],bucket->bucket_ring,PDEBUG);
92    }
93  #endif
94  pFalseReturn(p_Test(bucket->buckets[i], bucket->bucket_ring));
95  if (bucket->buckets_length[i] != pLength(bucket->buckets[i]))
96  {
97    dReportError("Bucket %d lengths difference should:%d has:%d",
98                 i, bucket->buckets_length[i], pLength(bucket->buckets[i]));
99  }
100  else if (i > 0 && (int) pLogLength(bucket->buckets_length[i]) > i)
101  {
102    dReportError("Bucket %d too long %d",
103                 i, bucket->buckets_length[i]);
104  }
105  if (i==0 && bucket->buckets_length[0] > 1)
106  {
107    dReportError("Bucket 0 too long");
108  }
109  return TRUE;
110}
111
112
113BOOLEAN kbTest(kBucket_pt bucket)
114{
115  #ifdef HAVE_COEF_BUCKETS
116  assume(bucket->coef[0]==NULL);
117  #endif
118  int i;
119  poly lm = bucket->buckets[0];
120
121  omCheckAddrBin(bucket, kBucket_bin);
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  //if (TEST_OPT_PROT) { Print("C(%d)",pl); }
432  kbTest(bucket);
433  return i;
434}
435
436void kBucketClear(kBucket_pt bucket, poly *p, int *length)
437{
438  int i = kBucketCanonicalize(bucket);
439  if (i > 0)
440  {
441  #ifdef USE_COEF_BUCKETS
442    MULTIPLY_BUCKET(bucket,i);
443    //bucket->coef[i]=NULL;
444  #endif
445    *p = bucket->buckets[i];
446    *length = bucket->buckets_length[i];
447    bucket->buckets[i] = NULL;
448    bucket->buckets_length[i] = 0;
449    bucket->buckets_used = 0;
450
451  }
452  else
453  {
454    *p = NULL;
455    *length = 0;
456  }
457}
458
459void kBucketSetLm(kBucket_pt bucket, poly lm)
460{
461  kBucketMergeLm(bucket);
462  pNext(lm) = NULL;
463  bucket->buckets[0] = lm;
464  bucket->buckets_length[0] = 1;
465}
466
467#else // HAVE_PSEUDO_BUCKETS
468
469void kBucketInit(kBucket_pt bucket, poly lm, int length)
470{
471  int i;
472
473  assume(bucket != NULL);
474  assume(length <= 0 || length == pLength(lm));
475
476  bucket->p = lm;
477  if (length <= 0) bucket->l = pLength(lm);
478  else bucket->l = length;
479
480}
481
482const poly kBucketGetLm(kBucket_pt bucket)
483{
484  return bucket->p;
485}
486
487poly kBucketExtractLm(kBucket_pt bucket)
488{
489  poly lm = bucket->p;
490  assume(pLength(bucket->p) == bucket->l);
491  pIter(bucket->p);
492  (bucket->l)--;
493  pNext(lm) = NULL;
494  return lm;
495}
496
497void kBucketClear(kBucket_pt bucket, poly *p, int *length)
498{
499  assume(pLength(bucket->p) == bucket->l);
500  *p = bucket->p;
501  *length = bucket->l;
502  bucket->p = NULL;
503  bucket->l = 0;
504}
505
506#endif // ! HAVE_PSEUDO_BUCKETS
507//////////////////////////////////////////////////////////////////////////
508///
509/// For changing the ring of the Bpoly to new_tailBin
510///
511void kBucketShallowCopyDelete(kBucket_pt bucket,
512                              ring new_tailRing, omBin new_tailBin,
513                              pShallowCopyDeleteProc p_shallow_copy_delete)
514{
515#ifndef HAVE_PSEUDO_BUCKETS
516  int i;
517
518  kBucketCanonicalize(bucket);
519  for (i=0; i<= bucket->buckets_used; i++)
520    if (bucket->buckets[i] != NULL)
521    {
522      MULTIPLY_BUCKET(bucket,i);
523      bucket->buckets[i] = p_shallow_copy_delete(bucket->buckets[i],
524                                                 bucket->bucket_ring,
525                                                 new_tailRing,
526                                                 new_tailBin);
527    }
528#else
529  bucket->p = p_shallow_copy_delete(p,
530                                    bucket_ring,
531                                    new_tailRing,
532                                    new_tailBin);
533#endif
534  bucket->bucket_ring = new_tailRing;
535}
536
537//////////////////////////////////////////////////////////////////////////
538///
539/// Bucket number i from bucket is out of length sync, resync
540///
541void kBucketAdjust(kBucket_pt bucket, int i) {
542
543  MULTIPLY_BUCKET(bucket,i);
544
545  int l1 = bucket->buckets_length[i];
546  poly p1 = bucket->buckets[i];
547  bucket->buckets[i] = NULL;
548  bucket->buckets_length[i] = 0;
549  i = pLogLength(l1);
550
551  while (bucket->buckets[i] != NULL)
552  {
553    //kbTest(bucket);
554    MULTIPLY_BUCKET(bucket,i);
555    p1 = p_Add_q(p1, bucket->buckets[i],
556                 l1, bucket->buckets_length[i], bucket->bucket_ring);
557    bucket->buckets[i] = NULL;
558    bucket->buckets_length[i] = 0;
559    i = pLogLength(l1);
560  }
561
562  bucket->buckets[i] = p1;
563  bucket->buckets_length[i]=l1;
564  if (i >= bucket->buckets_used)
565    bucket->buckets_used = i;
566  else
567    kBucketAdjustBucketsUsed(bucket);
568}
569
570//////////////////////////////////////////////////////////////////////////
571///
572/// Multiply Bucket by number ,i.e. Bpoly == n*Bpoly
573///
574void kBucket_Mult_n(kBucket_pt bucket, number n)
575{
576#ifndef HAVE_PSEUDO_BUCKETS
577  kbTest(bucket);
578  ring r=bucket->bucket_ring;
579  int i;
580
581  for (i=0; i<= bucket->buckets_used; i++)
582  {
583    if (bucket->buckets[i] != NULL)
584    {
585#ifdef USE_COEF_BUCKETS
586      if (i<coef_start)
587        bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
588#ifdef HAVE_RINGS
589        /* Frank Seelisch on March 11, 2010:
590           This looks a bit strange: The following "if" is indented
591           like the previous line of code. But coded as it is,
592           it should actually be two spaces less indented.
593           Question: Should the following "if" also only be
594           performed when "(i<coef_start)" is true?
595           For the time being, I leave it as it is. */
596        if (rField_is_Ring(r) && !(rField_is_Domain(r)))
597        {
598          bucket->buckets_length[i] = pLength(bucket->buckets[i]);
599          kBucketAdjust(bucket, i);
600        }
601#endif
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#ifdef HAVE_RINGS
614      if (rField_is_Ring(r) && !(rField_is_Domain(r)))
615      {
616        bucket->buckets_length[i] = pLength(bucket->buckets[i]);
617        kBucketAdjust(bucket, i);
618      }
619#endif
620#endif
621    }
622  }
623  kbTest(bucket);
624#else
625  bucket->p = p_Mult_nn(bucket->p, n, bucket->bucket_ring);
626#endif
627}
628
629
630//////////////////////////////////////////////////////////////////////////
631///
632/// Add to Bucket a poly ,i.e. Bpoly == q+Bpoly
633///
634void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
635{
636  if (q == NULL) return;
637  assume(*l <= 0 || pLength(q) == *l);
638
639  int i, l1;
640  ring r = bucket->bucket_ring;
641
642  if (*l <= 0)
643  {
644    l1 = pLength(q);
645    *l = l1;
646  }
647  else
648    l1 = *l;
649
650  kBucketMergeLm(bucket);
651  kbTest(bucket);
652  i = pLogLength(l1);
653
654  while (bucket->buckets[i] != NULL)
655  {
656    //MULTIPLY_BUCKET(bucket,i);
657  #ifdef USE_COEF_BUCKETS
658    if (bucket->coef[i]!=NULL)
659    {
660      q = p_Plus_mm_Mult_qq(q, bucket->coef[i], bucket->buckets[i],
661                 l1, bucket->buckets_length[i], r);
662      p_Delete(&bucket->coef[i],r);
663      p_Delete(&bucket->buckets[i],r);
664    }
665    else
666    q = p_Add_q(q, bucket->buckets[i],
667                 l1, bucket->buckets_length[i], r);
668  #else
669    q = p_Add_q(q, bucket->buckets[i],
670                 l1, bucket->buckets_length[i], r);
671  #endif
672    bucket->buckets[i] = NULL;
673    bucket->buckets_length[i] = 0;
674    i = pLogLength(l1);
675    assume(i<= MAX_BUCKET);
676    assume(bucket->buckets_used<= MAX_BUCKET);
677  }
678
679  kbTest(bucket);
680  bucket->buckets[i] = q;
681  bucket->buckets_length[i]=l1;
682  if (i >= bucket->buckets_used)
683    bucket->buckets_used = i;
684  else
685    kBucketAdjustBucketsUsed(bucket);
686  kbTest(bucket);
687}
688
689
690
691//////////////////////////////////////////////////////////////////////////
692///
693/// Bpoly == Bpoly - m*p; where m is a monom
694/// Does not destroy p and m
695/// assume (*l <= 0 || pLength(p) == *l)
696void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l,
697                            poly spNoether)
698{
699  assume(*l <= 0 || pLength(p) == *l);
700  int i, l1;
701  poly p1 = p;
702  poly last;
703  ring r = bucket->bucket_ring;
704
705  if (*l <= 0)
706  {
707    l1 = pLength(p1);
708    *l = l1;
709  }
710  else
711    l1 = *l;
712
713  if (m == NULL || p == NULL) return;
714
715#ifndef HAVE_PSEUDO_BUCKETS
716  kBucketMergeLm(bucket);
717  kbTest(bucket);
718  i = pLogLength(l1);
719
720  if ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
721  {
722    assume(pLength(bucket->buckets[i])==bucket->buckets_length[i]);
723//#ifdef USE_COEF_BUCKETS
724//     if(bucket->coef[i]!=NULL)
725//     {
726//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
727//       bucket->coef[i]=NULL;
728//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
729//                               bucket->buckets_length[i], l1,
730//                             spNoether, r);
731//     }
732//     else
733//#endif
734    MULTIPLY_BUCKET(bucket,i);
735    p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
736                            bucket->buckets_length[i], l1,
737                            spNoether, r);
738    l1 = bucket->buckets_length[i];
739    bucket->buckets[i] = NULL;
740    bucket->buckets_length[i] = 0;
741#ifdef HAVE_RINGS
742    if (rField_is_Ring(r) && !(rField_is_Domain(r)))
743    {
744      l1 = pLength(p1);
745      assume(pLength(p1) == l1);
746    }
747#endif
748    i = pLogLength(l1);
749  }
750  else
751  {
752    pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
753    if (spNoether != NULL)
754    {
755      l1 = -1;
756      p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r, last);
757      i = pLogLength(l1);
758    }
759    else {
760      p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
761#ifdef HAVE_RINGS
762      if (rField_is_Ring(r) && !(rField_is_Domain(r)))
763      {
764        l1 = pLength(p1);
765        i = pLogLength(l1);
766      }
767#endif
768    }
769    pSetCoeff0(m, n_Neg(pGetCoeff(m),r->cf));
770  }
771
772  while (bucket->buckets[i] != NULL)
773  {
774    //kbTest(bucket);
775    MULTIPLY_BUCKET(bucket,i);
776    p1 = p_Add_q(p1, bucket->buckets[i],
777                 l1, bucket->buckets_length[i], r);
778    bucket->buckets[i] = NULL;
779    bucket->buckets_length[i] = 0;
780    i = pLogLength(l1);
781  }
782
783  bucket->buckets[i] = p1;
784  bucket->buckets_length[i]=l1;
785  if (i >= bucket->buckets_used)
786    bucket->buckets_used = i;
787  else
788    kBucketAdjustBucketsUsed(bucket);
789#else // HAVE_PSEUDO_BUCKETS
790  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
791                               bucket->l, l1,
792                               spNoether, r);
793#endif
794}
795
796//////////////////////////////////////////////////////////////////////////
797///
798/// Bpoly == Bpoly + m*p; where m is a monom
799/// Does not destroy p and m
800/// assume (l <= 0 || pLength(p) == l)
801void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
802{
803    assume((!rIsPluralRing(bucket->bucket_ring))||p_IsConstant(m, bucket->bucket_ring));
804  assume(l <= 0 || pLength(p) == l);
805  int i, l1;
806  poly p1 = p;
807  poly last;
808  ring r = bucket->bucket_ring;
809
810  if (m == NULL || p == NULL) return;
811
812  if (l <= 0)
813  {
814    l1 = pLength(p1);
815    l = l1;
816  }
817  else
818    l1 = l;
819
820  kBucketMergeLm(bucket);
821  kbTest(bucket);
822  i = pLogLength(l1);
823  #ifdef USE_COEF_BUCKETS
824  number n=n_Init(1,r->cf);
825  #endif
826  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
827  {
828    //if (FALSE){
829    #ifdef USE_COEF_BUCKETS
830    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
831    {
832      number orig_coef=p_GetCoeff(bucket->coef[i],r);
833      //we take ownership:
834      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
835      number add_coef=n_Copy(p_GetCoeff(m,r),r);
836      number gcd=n_Gcd(add_coef, orig_coef,r);
837
838      if (!(n_IsOne(gcd,r)))
839      {
840        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
841        number add_coef2=n_IntDiv(add_coef, gcd,r);
842        n_Delete(&orig_coef,r);
843        n_Delete(&add_coef,r);
844        orig_coef=orig_coef2;
845        add_coef=add_coef2;
846
847        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
848        n_Delete(&n,r);
849        n=gcd;
850      }
851
852      //assume(n_IsOne(n,r));
853      number backup=p_GetCoeff(m,r);
854
855      p_SetCoeff0(m,add_coef,r);
856      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
857
858      n_Delete(&orig_coef,r);
859      p_Delete(&bucket->coef[i],r);
860
861      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
862                           bucket->buckets_length[i], l1, r);
863      l1=bucket->buckets_length[i];
864      bucket->buckets[i]=NULL;
865      bucket->buckets_length[i] = 0;
866      i = pLogLength(l1);
867      assume(l1==pLength(p1));
868
869      p_SetCoeff(m,backup,r); //deletes add_coef
870    }
871    else
872    #endif
873    {
874      MULTIPLY_BUCKET(bucket,i);
875      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
876                           bucket->buckets_length[i], l1, r);
877      l1 = bucket->buckets_length[i];
878      bucket->buckets[i] = NULL;
879      bucket->buckets_length[i] = 0;
880      i = pLogLength(l1);
881    }
882  }
883
884  else
885  {
886    #ifdef USE_COEF_BUCKETS
887    number swap_n=p_GetCoeff(m,r);
888
889    assume(n_IsOne(n,r));
890    p_SetCoeff0(m,n,r);
891    n=swap_n;
892    //p_SetCoeff0(n, swap_n, r);
893    //p_GetCoeff0(n, swap_n,r);
894    #endif
895    p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
896    #ifdef USE_COEF_BUCKETS
897    //m may not be changed
898    p_SetCoeff(m,n_Copy(n,r),r);
899    #endif
900  }
901
902
903  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
904  {
905    assume(i!=0);
906    #ifdef USE_COEF_BUCKETS
907    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
908    {
909      number orig_coef=p_GetCoeff(bucket->coef[i],r);
910      //we take ownership:
911      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
912      number add_coef=n_Copy(n,r);
913      number gcd=n_Gcd(add_coef, orig_coef,r);
914
915      if (!(n_IsOne(gcd,r)))
916      {
917        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
918        number add_coef2=n_IntDiv(add_coef, gcd,r);
919        n_Delete(&orig_coef,r);
920        n_Delete(&n,r);
921        n_Delete(&add_coef,r);
922        orig_coef=orig_coef2;
923        add_coef=add_coef2;
924        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
925        n=gcd;
926      }
927      //assume(n_IsOne(n,r));
928      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
929      p1=p_Mult_nn(p1,add_coef,r);
930
931      p1 = p_Add_q(p1, bucket->buckets[i],r);
932      l1=pLength(p1);
933
934      bucket->buckets[i]=NULL;
935      n_Delete(&orig_coef,r);
936      p_Delete(&bucket->coef[i],r);
937      //l1=bucket->buckets_length[i];
938      assume(l1==pLength(p1));
939    }
940    else
941    #endif
942    {
943      //don't do that, pull out gcd
944      #ifdef USE_COEF_BUCKETS
945      if(!(n_IsOne(n,r)))
946      {
947        p1=p_Mult_nn(p1, n, r);
948        n_Delete(&n,r);
949        n=n_Init(1,r);
950      }
951      #endif
952      MULTIPLY_BUCKET(bucket,i);
953      p1 = p_Add_q(p1, bucket->buckets[i],
954                 l1, bucket->buckets_length[i], r);
955      bucket->buckets[i] = NULL;
956      bucket->buckets_length[i] = 0;
957    }
958    i = pLogLength(l1);
959  }
960
961  bucket->buckets[i] = p1;
962#ifdef USE_COEF_BUCKETS
963  assume(bucket->coef[i]==NULL);
964
965  if (!(n_IsOne(n,r)))
966  {
967    bucket->coef[i]=p_NSet(n,r);
968  }
969  else
970  {
971    bucket->coef[i]=NULL;
972    n_Delete(&n,r);
973  }
974
975  if ((p1==NULL) && (bucket->coef[i]!=NULL))
976    p_Delete(&bucket->coef[i],r);
977#endif
978  bucket->buckets_length[i]=l1;
979  if (i >= bucket->buckets_used)
980    bucket->buckets_used = i;
981  else
982    kBucketAdjustBucketsUsed(bucket);
983
984  kbTest(bucket);
985}
986
987poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
988{
989  if (q == NULL) return append;
990  poly lm;
991  loop
992  {
993    lm = kBucketGetLm(bucket);
994    if (lm == NULL) return append;
995    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
996    {
997      lm = kBucketExtractLm(bucket);
998      pNext(append) = lm;
999      pIter(append);
1000    }
1001    else
1002    {
1003      return append;
1004    }
1005  }
1006}
1007
1008/////////////////////////////////////////////////////////////////////////////
1009//
1010// Extract all monomials from bucket with component comp
1011// Return as a polynomial *p with length *l
1012// In other words, afterwards
1013// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
1014// and components of monomials of *p are all 0
1015//
1016
1017// Hmm... for now I'm too lazy to implement those independent of currRing
1018// But better declare it extern than including polys.h
1019extern void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r);
1020
1021void kBucketTakeOutComp(kBucket_pt bucket,
1022                        long comp,
1023                        poly *r_p, int *l)
1024{
1025  poly p = NULL, q;
1026  int i, lp = 0, lq;
1027
1028#ifndef HAVE_PSEUDO_BUCKETS
1029  kBucketMergeLm(bucket);
1030  for (i=1; i<=bucket->buckets_used; i++)
1031  {
1032    if (bucket->buckets[i] != NULL)
1033    {
1034      MULTIPLY_BUCKET(bucket,i);
1035      p_TakeOutComp(&(bucket->buckets[i]), comp, &q, &lq, bucket->bucket_ring);
1036      if (q != NULL)
1037      {
1038        assume(pLength(q) == lq);
1039        bucket->buckets_length[i] -= lq;
1040        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
1041        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1042      }
1043    }
1044  }
1045  kBucketAdjustBucketsUsed(bucket);
1046#else
1047  p_TakeOutComp(&(bucket->p), comp, &p, &lp,bucket->bucket_ring);
1048  (bucket->l) -= lp;
1049#endif
1050  *r_p = p;
1051  *l = lp;
1052
1053  kbTest(bucket);
1054}
1055
1056/////////////////////////////////////////////////////////////////////////////
1057// Reduction of Bpoly with a given poly
1058//
1059
1060extern int ksCheckCoeff(number *a, number *b);
1061
1062number kBucketPolyRed(kBucket_pt bucket,
1063                      poly p1, int l1,
1064                      poly spNoether)
1065{
1066  ring r=bucket->bucket_ring;
1067  assume((!rIsPluralRing(r))||p_LmEqual(p1,kBucketGetLm(bucket), r));
1068  assume(p1 != NULL &&
1069         p_DivisibleBy(p1,  kBucketGetLm(bucket), r));
1070  assume(pLength(p1) == (int) l1);
1071
1072  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1073  BOOLEAN reset_vec=FALSE;
1074  number rn;
1075
1076  /* we shall reduce bucket=bn*lm+... by p1=an*t+a1 where t=lm(p1)
1077     and an,bn shall be defined further down only if lc(p1)!=1
1078     we already know: an|bn and t|lm */
1079  if(a1==NULL)
1080  {
1081    p_LmDelete(&lm, r);
1082    return n_Init(1,r->cf);
1083  }
1084
1085  if (! n_IsOne(pGetCoeff(p1),r->cf))
1086  {
1087    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1088//StringSetS("##### an = "); nWrite(an); PrintS(StringAppend("\n"));
1089//StringSetS("##### bn = "); nWrite(bn); PrintS(StringAppend("\n"));
1090    /* ksCheckCoeff: divide out gcd from an and bn: */
1091    int ct = ksCheckCoeff(&an, &bn,r->cf);
1092    /* the previous command returns ct=0 or ct=2 iff an!=1
1093       note: an is now 1 or -1 */
1094
1095    /* setup factor for p1 which cancels leading terms */
1096    p_SetCoeff(lm, bn, r);
1097    if ((ct == 0) || (ct == 2))
1098    {
1099      /* next line used to be here before but is WRONG:
1100      kBucket_Mult_n(bucket, an);
1101        its use would result in a wrong sign for the tail of bucket
1102        in the reduction */
1103
1104      /* correct factor for cancelation by changing sign if an=-1 */
1105      if (rField_is_Ring(r))
1106        lm = p_Mult_nn(lm, an, r);
1107      else
1108        kBucket_Mult_n(bucket, an);
1109    }
1110    rn = an;
1111  }
1112  else
1113  {
1114    rn = n_Init(1,r->cf);
1115  }
1116
1117  if (p_GetComp(p1, r) != p_GetComp(lm, r))
1118  {
1119    p_SetCompP(a1, p_GetComp(lm, r), r);
1120    reset_vec = TRUE;
1121    p_SetComp(lm, p_GetComp(p1, r), r);
1122    p_Setm(lm, r);
1123  }
1124
1125  p_ExpVectorSub(lm, p1, r);
1126  l1--;
1127
1128  assume(l1==pLength(a1));
1129  BOOLEAN backuped=FALSE;
1130  number coef;
1131  #if 0
1132  //@Viktor, don't ignore coefficients on monomials
1133  if(l1==1) {
1134
1135    //if (rField_is_Q(r)) {
1136      //avoid this for function fields, as gcds are expensive at the moment
1137
1138
1139      coef=p_GetCoeff(a1,r);
1140      lm=p_Mult_nn(lm, coef, r);
1141      p_SetCoeff0(a1, n_Init(1,r), r);
1142      backuped=TRUE;
1143      //WARNING: not thread_safe
1144    //deletes coef as side effect
1145    //}
1146  }
1147  #endif
1148
1149  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1150
1151  if (backuped)
1152    p_SetCoeff0(a1,coef,r);
1153  p_LmDelete(&lm, r);
1154  if (reset_vec) p_SetCompP(a1, 0, r);
1155  kbTest(bucket);
1156  return rn;
1157}
1158static BOOLEAN nIsPseudoUnit(number n, ring r)
1159{
1160  if (rField_is_Zp(r))
1161    return TRUE;
1162
1163  if (rParameter(r)==NULL)
1164  {
1165    return (n_Size(n,r->cf)==1);
1166  }
1167  //if (r->parameter!=NULL)
1168  return (n_IsOne(n,r->cf) || n_IsMOne(n,r->cf));
1169}
1170
1171#ifndef USE_COEF_BUCKETS
1172void kBucketSimpleContent(kBucket_pt) {}
1173#else
1174void kBucketSimpleContent(kBucket_pt bucket)
1175{
1176  ring r=bucket->bucket_ring;
1177  int i;
1178  //PrintS("HHHHHHHHHHHHH");
1179  for (i=0;i<=MAX_BUCKET;i++)
1180  {
1181    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1182    //    PrintS("H2H2H2");
1183    if (i==0)
1184    {
1185      assume(bucket->buckets[i]==NULL);
1186    }
1187    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1188      return;
1189  }
1190  for (i=0;i<=MAX_BUCKET;i++)
1191  {
1192    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1193    //    PrintS("H2H2H2");
1194    if (i==0)
1195    {
1196      assume(bucket->buckets[i]==NULL);
1197    }
1198    if ((bucket->buckets[i]!=NULL)
1199    && (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1200      return;
1201  }
1202  //return;
1203
1204  number coef=n_Init(0,r);
1205  //ATTENTION: will not work correct for GB over ring
1206  //if (TEST_OPT_PROT)
1207  //    PrintS("CCCCCCCCCCCCC");
1208  for (i=MAX_BUCKET;i>=0;i--)
1209  {
1210    if (i==0)
1211    {
1212      assume(bucket->buckets[i]==NULL);
1213    }
1214    if (bucket->buckets[i]!=NULL)
1215    {
1216      assume(bucket->coef[i]!=NULL);
1217      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1218
1219      //in this way it should crash on programming errors, yeah
1220      number temp=n_Gcd(coef, pGetCoeff(bucket->coef[i]),r);
1221      n_Delete(&coef,r );
1222      coef=temp;
1223      if (nIsPseudoUnit(coef,r))
1224      {
1225        n_Delete(&coef,r);
1226        return;
1227      }
1228      assume(!(n_IsZero(coef,r)));
1229    }
1230  }
1231  if (n_IsZero(coef,r))
1232  {
1233    n_Delete(&coef,r);
1234    return;
1235  }
1236  if (TEST_OPT_PROT)
1237    PrintS("S");
1238  for(i=0;i<=MAX_BUCKET;i++)
1239  {
1240    if (bucket->buckets[i]!=NULL)
1241    {
1242      assume(!(n_IsZero(coef,r)));
1243      assume(bucket->coef[i]!=NULL);
1244      number lc=p_GetCoeff(bucket->coef[i],r);
1245      p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1246      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1247    }
1248  }
1249  n_Delete(&coef,r);
1250}
1251#endif
1252
1253
1254poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1255{
1256  assume(bucket->buckets[i]!=NULL);
1257
1258  poly p=bucket->buckets[i];
1259  bucket->buckets_length[i]--;
1260#ifdef USE_COEF_BUCKETS
1261  ring r=bucket->bucket_ring;
1262  if (bucket->coef[i]!=NULL)
1263  {
1264    poly next=pNext(p);
1265    if (next==NULL)
1266    {
1267      MULTIPLY_BUCKET(bucket,i);
1268      p=bucket->buckets[i];
1269      bucket->buckets[i]=NULL;
1270      return p;
1271    }
1272    else
1273    {
1274      bucket->buckets[i]=next;
1275      number c=p_GetCoeff(bucket->coef[i],r);
1276      pNext(p)=NULL;
1277      p=p_Mult_nn(p,c,r);
1278      assume(p!=NULL);
1279      return p;
1280    }
1281  }
1282  else
1283#endif
1284  {
1285    bucket->buckets[i]=pNext(bucket->buckets[i]);
1286    pNext(p)=NULL;
1287    assume(p!=NULL);
1288    return p;
1289  }
1290}
1291
1292/*
1293* input - output: a, b
1294* returns:
1295*   a := a/gcd(a,b), b := b/gcd(a,b)
1296*   and return value
1297*       0  ->  a != 1,  b != 1
1298*       1  ->  a == 1,  b != 1
1299*       2  ->  a != 1,  b == 1
1300*       3  ->  a == 1,  b == 1
1301*   this value is used to control the spolys
1302*/
1303int ksCheckCoeff(number *a, number *b, const coeffs r)
1304{
1305  int c = 0;
1306  number an = *a, bn = *b;
1307  n_Test(an,r);
1308  n_Test(bn,r);
1309
1310  number cn = n_Gcd(an, bn, r);
1311
1312  if(n_IsOne(cn, r))
1313  {
1314    an = n_Copy(an, r);
1315    bn = n_Copy(bn, r);
1316  }
1317  else
1318  {
1319    an = n_IntDiv(an, cn, r);
1320    bn = n_IntDiv(bn, cn, r);
1321  }
1322  n_Delete(&cn, r);
1323  if (n_IsOne(an, r))
1324  {
1325    c = 1;
1326  }
1327  if (n_IsOne(bn, r))
1328  {
1329    c += 2;
1330  }
1331  *a = an;
1332  *b = bn;
1333  return c;
1334}
1335
Note: See TracBrowser for help on using the repository browser.