source: git/kernel/tgb_internal.h @ 44fac1

jengelh-datetimespielwiese
Last change on this file since 44fac1 was 44fac1, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: omalloc -> omAlloc from master
  • Property mode set to 100644
File size: 49.4 KB
Line 
1#ifndef TGB_INTERNAL_H
2#define TGB_INTERNAL_H
3//!\file tgb_internal.h
4/****************************************
5*  Computer Algebra System SINGULAR     *
6****************************************/
7/*
8 * ABSTRACT: tgb internal .h file
9*/
10#define USE_NORO 1
11
12#include <omalloc/omalloc.h>
13
14//#define TGB_DEBUG
15#define FULLREDUCTIONS
16#define HANS_IDEA
17//#define HALFREDUCTIONS
18//#define HEAD_BIN
19//#define HOMOGENEOUS_EXAMPLE
20#define REDTAIL_S
21#define PAR_N 100
22#define PAR_N_F4 5000
23#define AC_NEW_MIN 2
24#define AC_FLATTEN 1
25
26//#define FIND_DETERMINISTIC
27//#define REDTAIL_PROT
28//#define QUICK_SPOLY_TEST
29#ifdef USE_NORO
30#define NORO_CACHE 1
31#define NORO_SPARSE_ROWS_PRE 1
32#define NORO_NON_POLY 1
33#include <algorithm>
34#endif
35#ifdef NORO_CACHE
36//#include <map>
37#include <vector>
38#endif
39#ifdef HAVE_BOOST_DYNAMIC_BITSET_HPP
40#define  HAVE_BOOST 1
41#endif
42//#define HAVE_BOOST 1
43//#define USE_STDVECBOOL 1
44#ifdef HAVE_BOOST
45#include <vector>
46using boost::dynamic_bitset;
47using std::vector;
48#endif
49#ifdef USE_STDVECBOOL
50#include <vector>
51using std::vector;
52#endif
53#include <stdlib.h>
54
55#include <polys/monomials/p_polys.h>
56
57#include <kernel/ideals.h>
58#include <polys/monomials/ring.h>
59#include <kernel/febase.h>
60#include <misc/options.h>
61#include <kernel/polys.h>
62#include <kernel/modulop.h>
63#include <kernel/kutil.h>
64#include <kernel/kInline.h>
65#include <kernel/kstd1.h>
66#include <polys/kbuckets.h>
67
68class PolySimple
69{
70public:
71  PolySimple(poly p)
72  {
73    impl=p;
74  }
75  PolySimple()
76  {
77    impl=NULL;
78  }
79  PolySimple(const PolySimple& a)
80  {
81    //impl=p_Copy(a.impl,currRing);
82    impl=a.impl;
83  }
84  PolySimple& operator=(const PolySimple& p2)
85  {
86    //p_Delete(&impl,currRing);
87    //impl=p_Copy(p2.impl,currRing);
88    impl=p2.impl;
89    return *this;
90  }
91  ~PolySimple()
92  {
93    //p_Delete(&impl,currRing);
94  }
95  bool operator< (const PolySimple& other) const
96  {
97    return pLmCmp(impl,other.impl)<0;
98  }
99  bool operator==(const PolySimple& other)
100  {
101    return pLmEqual(impl,other.impl);
102  }
103  poly impl;
104
105};
106template<class number_type> class DataNoroCacheNode;
107/*class MonRedRes{
108public:
109  poly p;
110  number coef;
111  BOOLEAN changed;
112  int len;
113  BOOLEAN onlyBorrowed;
114  bool operator<(const MonRedRes& other) const{
115    int cmp=p_LmCmp(p,other.p,currRing);
116    if ((cmp<0)||((cmp==0)&&((onlyBorrowed)&&(!(other.onlyBorrowed))))){
117      return true;
118    } else return false;
119  }
120  DataNoroCacheNode* ref;
121  MonRedRes(){
122    ref=NULL;
123    p=NULL;
124  }
125};*/
126template <class number_type> class MonRedResNP
127{
128public:
129  number coef;
130
131
132  DataNoroCacheNode<number_type>* ref;
133  MonRedResNP()
134  {
135    ref=NULL;
136  }
137};
138struct sorted_pair_node
139{
140  //criterium, which is stable 0. small lcm 1. small i 2. small j
141  wlen_type expected_length;
142  poly lcm_of_lm;
143  int i;
144  int j;
145  int deg;
146
147
148};
149#ifdef NORO_CACHE
150#ifndef NORO_NON_POLY
151class NoroPlaceHolder
152{
153public:
154  DataNoroCacheNode* ref;
155  number coef;
156};
157#endif
158#endif
159//static ideal debug_Ideal;
160
161
162struct poly_list_node
163{
164  poly p;
165  poly_list_node* next;
166};
167
168struct int_pair_node
169{
170  int_pair_node* next;
171  int a;
172  int b;
173};
174struct monom_poly
175{
176  poly m;
177  poly f;
178};
179struct mp_array_list
180{
181  monom_poly* mp;
182  int size;
183  mp_array_list* next;
184};
185
186
187struct poly_array_list
188{
189  poly* p;
190  int size;
191  poly_array_list* next;
192};
193class slimgb_alg
194{
195  public:
196    slimgb_alg(ideal I, int syz_comp,BOOLEAN F4,int deg_pos);
197                void introduceDelayedPairs(poly* pa,int s);
198    virtual ~slimgb_alg();
199    void cleanDegs(int lower, int upper);
200#ifndef HAVE_BOOST
201#ifdef USE_STDVECBOOL
202  vector<vector<bool> > states;
203#else
204  char** states;
205#endif
206#else
207  vector<dynamic_bitset<> > states;
208#endif
209  ideal add_later;
210  ideal S;
211  ring r;
212  int* lengths;
213  wlen_type* weighted_lengths;
214  long* short_Exps;
215  kStrategy strat;
216  int* T_deg;
217  int* T_deg_full;
218  poly tmp_lm;
219  poly* tmp_pair_lm;
220  sorted_pair_node** tmp_spn;
221  poly* expandS;
222  poly* gcd_of_terms;
223  int_pair_node* soon_free;
224  sorted_pair_node** apairs;
225  #if 0
226  BOOLEAN* modifiedS;
227  #endif
228  #ifdef TGB_RESORT_PAIRS
229  bool* replaced;
230  #endif
231  poly_list_node* to_destroy;
232  //for F4
233  mp_array_list* F;
234  poly_array_list* F_minus;
235
236  //end for F4
237#ifdef HEAD_BIN
238  omBin   HeadBin;
239#endif
240  unsigned int reduction_steps;
241  int n;
242  //! array_lengths should be greater equal n;
243  int syz_comp;
244  int array_lengths;
245  int normal_forms;
246  int current_degree;
247  int Rcounter;
248  int last_index;
249  int max_pairs;
250  int pair_top;
251  int easy_product_crit;
252  int extended_product_crit;
253  int average_length;
254  int lastDpBlockStart;
255  int lastCleanedDeg;
256  int deg_pos;
257  BOOLEAN use_noro;
258  BOOLEAN use_noro_last_block;
259  BOOLEAN isDifficultField;
260  BOOLEAN completed;
261  BOOLEAN is_homog;
262  BOOLEAN tailReductions;
263  BOOLEAN eliminationProblem;
264  BOOLEAN F4_mode;
265  BOOLEAN nc;
266  #ifdef TGB_RESORT_PAIRS
267  BOOLEAN used_b;
268  #endif
269  unsigned long pTotaldegree(poly p)
270  {
271      pTest(p);
272      //assume(pDeg(p,r)==::p_Totaldegree(p,r));
273      assume(((unsigned long)::p_Totaldegree(p,r))==p->exp[deg_pos]);
274      return p->exp[deg_pos];
275      //return ::pTotaldegree(p,this->r);
276  }
277  int pTotaldegree_full(poly p)
278  {
279    int rr=0;
280    while(p)
281    {
282      int d=this->pTotaldegree(p);
283      rr=si_max(rr,d);
284      pIter(p);
285    }
286    return rr;
287  }
288};
289class red_object
290{
291 public:
292  kBucket_pt bucket;
293  poly p;
294  unsigned long sev;
295  int sugar;
296  void flatten();
297  void validate();
298  wlen_type initial_quality;
299  void adjust_coefs(number c_r, number c_ac_r);
300  wlen_type guess_quality(slimgb_alg* c);
301  int clear_to_poly();
302  void canonicalize();
303};
304
305
306enum calc_state
307  {
308    UNCALCULATED,
309    HASTREP//,
310    //UNIMPORTANT,
311    //SOONTREP
312  };
313//static BOOLEAN pair_cmp(sorted_pair_node* a,sorted_pair_node* b);
314template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set);
315static int add_to_reductors(slimgb_alg* c, poly h, int len, int ecart, BOOLEAN simplified=FALSE);
316static int bucket_guess(kBucket* bucket);
317static poly redNFTail (poly h,const int sl,kStrategy strat, int len);
318static poly redNF2 (poly h,slimgb_alg* c , int &len, number&  m,int n=0);
319void free_sorted_pair_node(sorted_pair_node* s, ring r);
320static void shorten_tails(slimgb_alg* c, poly monom);
321static void replace_pair(int & i, int & j, slimgb_alg* c);
322//static sorted_pair_node** add_to_basis(poly h, int i, int j,slimgb_alg* c, int* ip=NULL);
323//static void do_this_spoly_stuff(int i,int j,slimgb_alg* c);
324//ideal t_rep_gb(ring r,ideal arg_I);
325ideal do_t_rep_gb(ring r,ideal arg_I, int syz_comp, BOOLEAN F4_mode,int deg_pos);
326static BOOLEAN has_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* state);
327//static int* make_connections(int from, poly bound, slimgb_alg* c);
328static int* make_connections(int from, int to, poly bound, slimgb_alg* c);
329void now_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
330//static void soon_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
331//static int pLcmDeg(poly a, poly b);
332static int simple_posInS (kStrategy strat, poly p,int len, wlen_type wlen);
333//static BOOLEAN find_next_pair(slimgb_alg* c, BOOLEAN go_higher=TRUE);
334
335//static sorted_pair_node* pop_pair(slimgb_alg* c);
336//static BOOLEAN no_pairs(slimgb_alg* c);
337void clean_top_of_pair_list(slimgb_alg* c);
338static void super_clean_top_of_pair_list(slimgb_alg* c);
339static BOOLEAN state_is(calc_state state, const int & i, const int & j, slimgb_alg* c);
340static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, slimgb_alg* c=NULL);
341static int tgb_pair_better_gen(const void* ap,const void* bp);
342static poly redTailShort(poly h, kStrategy strat);
343static poly gcd_of_terms(poly p, ring r);
344static BOOLEAN extended_product_criterion(poly p1, poly gcd1, poly p2, poly gcd2, slimgb_alg* c);
345//static poly kBucketGcd(kBucket* b, ring r);
346static void multi_reduction(red_object* los, int & losl, slimgb_alg* c);
347int slim_nsize(number n, ring r);
348sorted_pair_node* quick_pop_pair(slimgb_alg* c);
349sorted_pair_node* top_pair(slimgb_alg* c);
350sorted_pair_node** add_to_basis_ideal_quotient(poly h, slimgb_alg* c, int* ip);//, BOOLEAN new_pairs=TRUE);
351sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,slimgb_alg* c);
352int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj);
353int tgb_pair_better_gen2(const void* ap,const void* bp);
354int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev);
355//static int quality(poly p, int len, slimgb_alg* c);
356/**
357   makes on each red_object in a region a single_step
358 **/
359class reduction_step
360{
361 public:
362  /// we assume hat all occuring red_objects have same lm, and all
363  /// occ. lm's in r[l...u] are the same, only reductor does not occur
364  virtual void reduce(red_object* r, int l, int u);
365  //int reduction_id;
366  virtual ~reduction_step();
367  slimgb_alg* c;
368  int reduction_id;
369};
370class simple_reducer:public reduction_step
371{
372 public:
373  poly p;
374  kBucket_pt fill_back;
375  int p_len;
376  int reducer_deg;
377  simple_reducer(poly pp, int pp_len,int pp_reducer_deg, slimgb_alg* pp_c =NULL)
378  {
379    this->p=pp;
380    this->reducer_deg=pp_reducer_deg;
381    assume(pp_len==pLength(pp));
382    this->p_len=pp_len;
383    this->c=pp_c;
384  }
385  virtual void pre_reduce(red_object* r, int l, int u);
386  virtual void reduce(red_object* r, int l, int u);
387  ~simple_reducer();
388
389
390  virtual void do_reduce(red_object & ro);
391};
392
393//class sum_canceling_reducer:public reduction_step {
394//  void reduce(red_object* r, int l, int u);
395//};
396struct find_erg
397{
398  poly expand;
399  int expand_length;
400  int to_reduce_u;
401  int to_reduce_l;
402  int reduce_by;//index of reductor
403  BOOLEAN fromS;//else from los
404
405};
406
407static void multi_reduce_step(find_erg & erg, red_object* r, slimgb_alg* c);
408//static void finalize_reduction_step(reduction_step* r);
409
410template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set)
411{
412  //Print("POSHELER:%d",sizeof(wlen_type));
413  int length=strat->sl;
414  int i;
415  int an = 0;
416  int en= length;
417
418  if ((len>setL[length])
419      || ((len==setL[length]) && (pLmCmp(set[length],p)== -1)))
420    return length+1;
421
422  loop
423  {
424    if (an >= en-1)
425    {
426      if ((len<setL[an])
427          || ((len==setL[an]) && (pLmCmp(set[an],p) == 1))) return an;
428      return en;
429    }
430    i=(an+en) / 2;
431    if ((len<setL[i])
432        || ((len==setL[i]) && (pLmCmp(set[i],p) == 1))) en=i;
433    //else if ((len>setL[i])
434    //|| ((len==setL[i]) && (pLmCmp(set[i],p) == -1))) an=i;
435    else an=i;
436  }
437
438}
439#ifdef NORO_CACHE
440#define slim_prec_cast(a) (unsigned int) (unsigned long) (a)
441#define F4mat_to_number_type(a) (number_type) slim_prec_cast(a)
442typedef unsigned short tgb_uint16;
443typedef unsigned char tgb_uint8;
444typedef unsigned int tgb_uint32;
445class NoroCacheNode
446{
447public:
448  NoroCacheNode** branches;
449  int branches_len;
450
451
452  NoroCacheNode()
453  {
454    branches=NULL;
455    branches_len=0;
456
457  }
458  NoroCacheNode* setNode(int branch, NoroCacheNode* node)
459  {
460    if (branch>=branches_len)
461    {
462      if (branches==NULL)
463      {
464        branches_len=branch+1;
465        branches_len=si_max(branches_len,3);
466        branches=(NoroCacheNode**) omAlloc(branches_len*sizeof(NoroCacheNode*));
467        int i;
468        for(i=0;i<branches_len;i++)
469        {
470          branches[i]=NULL;
471        }
472      }
473      else
474      {
475        int branches_len_old=branches_len;
476        branches_len=branch+1;
477        branches=(NoroCacheNode**) omrealloc(branches,branches_len*sizeof(NoroCacheNode*));
478        int i;
479        for(i=branches_len_old;i<branches_len;i++)
480        {
481          branches[i]=NULL;
482        }
483      }
484    }
485    assume(branches[branch]==NULL);
486    branches[branch]=node;
487    return node;
488  }
489  NoroCacheNode* getBranch(int branch)
490  {
491    if (branch<branches_len) return branches[branch];
492    return NULL;
493  }
494  virtual ~NoroCacheNode()
495  {
496    int i;
497    for(i=0;i<branches_len;i++)
498    {
499      delete branches[i];
500    }
501    omfree(branches);
502  }
503  NoroCacheNode* getOrInsertBranch(int branch)
504  {
505    if ((branch<branches_len)&&(branches[branch]))
506      return branches[branch];
507    else
508    {
509      return setNode(branch,new NoroCacheNode());
510    }
511  }
512};
513class DenseRow{
514public:
515  number* array;
516  int begin;
517  int end;
518  DenseRow()
519  {
520    array=NULL;
521  }
522  ~DenseRow()
523  {
524    omfree(array);
525  }
526};
527template <class number_type> class SparseRow
528{
529public:
530  int* idx_array;
531  number_type* coef_array;
532  int len;
533  SparseRow()
534  {
535    len=0;
536    idx_array=NULL;
537    coef_array=NULL;
538  }
539  SparseRow<number_type>(int n)
540  {
541    len=n;
542    idx_array=(int*) omAlloc(n*sizeof(int));
543    coef_array=(number_type*) omAlloc(n*sizeof(number_type));
544  }
545  SparseRow<number_type>(int n, const number_type* source)
546  {
547    len=n;
548    idx_array=NULL;
549    coef_array=(number_type*) omAlloc(n*sizeof(number_type));
550    memcpy(coef_array,source,n*sizeof(number_type));
551  }
552  ~SparseRow<number_type>()
553  {
554    omfree(idx_array);
555    omfree(coef_array);
556  }
557};
558
559template <class number_type> class DataNoroCacheNode:public NoroCacheNode
560{
561public:
562
563  int value_len;
564  poly value_poly;
565  #ifdef NORO_SPARSE_ROWS_PRE
566  SparseRow<number_type>* row;
567  #else
568  DenseRow* row;
569  #endif
570  int term_index;
571  DataNoroCacheNode(poly p, int len)
572  {
573    value_len=len;
574    value_poly=p;
575    row=NULL;
576    term_index=-1;
577  }
578  #ifdef NORO_SPARSE_ROWS_PRE
579  DataNoroCacheNode(SparseRow<number_type>* row)
580  {
581    if (row!=NULL)
582      value_len=row->len;
583    else
584      value_len=0;
585    value_poly=NULL;
586    this->row=row;
587    term_index=-1;
588  }
589  #endif
590  ~DataNoroCacheNode()
591  {
592    //p_Delete(&value_poly,currRing);
593    if (row) delete row;
594  }
595};
596template <class number_type> class TermNoroDataNode
597{
598public:
599  DataNoroCacheNode<number_type>* node;
600  poly t;
601};
602
603template <class number_type> class NoroCache
604{
605public:
606  poly temp_term;
607#ifndef NORO_NON_POLY
608  void evaluatePlaceHolder(number* row,std::vector<NoroPlaceHolder>& place_holders);
609  void evaluateRows();
610  void evaluateRows(int level, NoroCacheNode* node);
611#endif
612  void collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type>* >& res);
613  void collectIrreducibleMonomials(int level,  NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>* >& res);
614
615#ifdef NORO_RED_ARRAY_RESERVER
616  int reserved;
617  poly* recursionPolyBuffer;
618#endif
619  static const int backLinkCode=-222;
620  DataNoroCacheNode<number_type>* insert(poly term, poly nf, int len)
621  {
622    //assume(impl.find(p_Copy(term,currRing))==impl.end());
623    //assume(len==pLength(nf));
624    assume(n_IsOne(p_GetCoeff(term,currRing),currRing->cf));
625    if (term==nf)
626    {
627      term=p_Copy(term,currRing);
628
629      ressources.push_back(term);
630      nIrreducibleMonomials++;
631      return treeInsertBackLink(term);
632
633    }
634    else
635    {
636      if (nf)
637      {
638        //nf=p_Copy(nf,currRing);
639        assume(p_LmCmp(nf,term,currRing)==-1);
640        ressources.push_back(nf);
641      }
642      return treeInsert(term,nf,len);
643
644    }
645
646    //impl[term]=std::pair<PolySimple,int> (nf,len);
647  }
648  #ifdef NORO_SPARSE_ROWS_PRE
649  DataNoroCacheNode<number_type>* insert(poly term, SparseRow<number_type>* srow)
650  {
651    //assume(impl.find(p_Copy(term,currRing))==impl.end());
652    //assume(len==pLength(nf));
653
654      return treeInsert(term,srow);
655
656
657    //impl[term]=std::pair<PolySimple,int> (nf,len);
658  }
659  #endif
660  DataNoroCacheNode<number_type>* insertAndTransferOwnerShip(poly t, ring r)
661  {
662    ressources.push_back(t);
663    DataNoroCacheNode<number_type>* res=treeInsertBackLink(t);
664    res->term_index=nIrreducibleMonomials;
665    nIrreducibleMonomials++;
666    return res;
667  }
668  poly lookup(poly term, BOOLEAN& succ, int & len);
669  DataNoroCacheNode<number_type>* getCacheReference(poly term);
670  NoroCache()
671  {
672    buffer=NULL;
673#ifdef NORO_RED_ARRAY_RESERVER
674    reserved=0;
675    recursionPolyBuffer=(poly*)omAlloc(1000000*sizeof(poly));
676#endif
677    nIrreducibleMonomials=0;
678    nReducibleMonomials=0;
679    temp_term=pOne();
680    tempBufferSize=3000;
681    tempBuffer=omAlloc(tempBufferSize);
682  }
683  void ensureTempBufferSize(size_t size)
684  {
685    if (tempBufferSize<size)
686    {
687      tempBufferSize=2*size;
688      omFree(tempBuffer);
689      tempBuffer=omAlloc(tempBufferSize);
690    }
691  }
692#ifdef NORO_RED_ARRAY_RESERVER
693  poly* reserve(int n)
694  {
695    poly* res=recursionPolyBuffer+reserved;
696    reserved+=n;
697    return res;
698  }
699  void free(int n)
700  {
701    reserved-=n;
702  }
703#endif
704  ~NoroCache()
705  {
706    int s=ressources.size();
707    int i;
708    for(i=0;i<s;i++)
709    {
710      p_Delete(&ressources[i].impl,currRing);
711    }
712    p_Delete(&temp_term,currRing);
713#ifdef NORO_RED_ARRAY_RESERVER
714    omfree(recursionPolyBuffer);
715#endif
716   omFree(tempBuffer);
717  }
718
719  int nIrreducibleMonomials;
720  int nReducibleMonomials;
721  void* tempBuffer;
722  size_t tempBufferSize;
723protected:
724  DataNoroCacheNode<number_type>* treeInsert(poly term,poly nf,int len)
725  {
726    int i;
727    nReducibleMonomials++;
728    int nvars=(currRing->N);
729    NoroCacheNode* parent=&root;
730    for(i=1;i<nvars;i++)
731    {
732      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
733    }
734    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(nf,len));
735  }
736  #ifdef NORO_SPARSE_ROWS_PRE
737  DataNoroCacheNode<number_type>* treeInsert(poly term,SparseRow<number_type>* srow)
738  {
739    int i;
740    nReducibleMonomials++;
741    int nvars=(currRing->N);
742    NoroCacheNode* parent=&root;
743    for(i=1;i<nvars;i++)
744    {
745      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
746    }
747    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(srow));
748  }
749  #endif
750  DataNoroCacheNode<number_type>* treeInsertBackLink(poly term)
751  {
752    int i;
753    int nvars=(currRing->N);
754    NoroCacheNode* parent=&root;
755    for(i=1;i<nvars;i++)
756    {
757      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
758    }
759    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(term,backLinkCode));
760  }
761
762  //@TODO descruct nodes;
763  typedef std::vector<PolySimple> poly_vec;
764  poly_vec ressources;
765  //typedef std::map<PolySimple,std::pair<PolySimple,int> > cache_map;
766  //cache_map impl;
767  NoroCacheNode root;
768  number* buffer;
769};
770template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c);
771template<class number_type> MonRedResNP<number_type> noro_red_mon_to_non_poly(poly t,  NoroCache<number_type> * cache,slimgb_alg* c)
772{
773  MonRedResNP<number_type> res_holder;
774
775
776    DataNoroCacheNode<number_type>* ref=cache->getCacheReference(t);
777    if (ref!=NULL)
778    {
779      res_holder.coef=p_GetCoeff(t,c->r);
780
781      res_holder.ref=ref;
782      p_Delete(&t,c->r);
783      return res_holder;
784    }
785
786  unsigned long sev=p_GetShortExpVector(t,currRing);
787  int i=kFindDivisibleByInS_easy(c->strat,t,sev);
788  if (i>=0)
789  {
790    number coef_bak=p_GetCoeff(t,c->r);
791
792    p_SetCoeff(t,n_Init(1,c->r->cf),c->r);
793    assume(n_IsOne(p_GetCoeff(c->strat->S[i],c->r),c->r->cf));
794    number coefstrat=p_GetCoeff(c->strat->S[i],c->r);
795
796
797    poly exp_diff=cache->temp_term;
798    p_ExpVectorDiff(exp_diff,t,c->strat->S[i],c->r);
799    p_SetCoeff(exp_diff,n_Neg(n_Invers(coefstrat,c->r->cf),c->r->cf),c->r);
800    p_Setm(exp_diff,c->r);
801    assume(c->strat->S[i]!=NULL);
802
803    poly res;
804    res=pp_Mult_mm(pNext(c->strat->S[i]),exp_diff,c->r);
805
806    int len=c->strat->lenS[i]-1;
807    SparseRow<number_type>* srow;
808    srow=noro_red_to_non_poly_t<number_type>(res,len,cache,c);
809    ref=cache->insert(t,srow);
810    p_Delete(&t,c->r);
811
812
813    res_holder.coef=coef_bak;
814    res_holder.ref=ref;
815    return res_holder;
816
817  } else {
818    number coef_bak=p_GetCoeff(t,c->r);
819    number one=n_Init(1, c->r->cf);
820    p_SetCoeff(t,one,c->r);
821
822    res_holder.ref=cache->insertAndTransferOwnerShip(t,c->r);
823    assume(res_holder.ref!=NULL);
824    res_holder.coef=coef_bak;
825
826    return res_holder;
827
828  }
829
830}
831/*
832poly tree_add(poly* a,int begin, int end,ring r)
833{
834  int d=end-begin;
835  switch(d)
836  {
837    case 0:
838      return NULL;
839    case 1:
840      return a[begin];
841    case 2:
842      return p_Add_q(a[begin],a[begin+1],r);
843    default:
844      int s=d/2;
845      return p_Add_q(tree_add(a,begin,begin+s,r),tree_add(a,begin+s,end,r),r);
846  }
847}
848*/
849#ifdef __GNUC__
850#define LIKELY(expression) (__builtin_expect(!!(expression), 1))
851#define UNLIKELY(expression) (__builtin_expect(!!(expression), 0))
852#else
853#define LIKELY(expression) (expression)
854#define UNLIKELY(expression) (expression)
855#endif
856
857template<class number_type> SparseRow<number_type>* convert_to_sparse_row(number_type* temp_array,int temp_size,int non_zeros){
858SparseRow<number_type>* res=new SparseRow<number_type>(non_zeros);
859//int pos=0;
860//Print("denseness:%f\n",((double) non_zeros/(double) temp_size));
861number_type* it_coef=res->coef_array;
862int* it_idx=res->idx_array;
863#if 0
864for(i=0;i<cache->nIrreducibleMonomials;i++){
865  if (!(0==temp_array[i])){
866
867    res->idx_array[pos]=i;
868    res->coef_array[pos]=temp_array[i];
869
870    pos++;
871    non_zeros--;
872    if (non_zeros==0) break;
873  }
874
875}
876#else
877int64* start=(int64*) ((void*)temp_array);
878int64* end;
879const int multiple=sizeof(int64)/sizeof(number_type);
880if (temp_size==0) end=start;
881
882else
883{
884  int temp_size_rounded=temp_size+(multiple-(temp_size%multiple));
885  assume(temp_size_rounded>=temp_size);
886  assume(temp_size_rounded%multiple==0);
887  assume(temp_size_rounded<temp_size+multiple);
888  number_type* nt_end=temp_array+temp_size_rounded;
889  end=(int64*)((void*)nt_end);
890}
891int64* it=start;
892while(it!=end)
893{
894  if UNLIKELY((*it)!=0)
895  {
896    int small_i;
897    const int temp_index=((number_type*)((void*) it))-temp_array;
898    const int bound=temp_index+multiple;
899    number_type c;
900    for(small_i=temp_index;small_i<bound;small_i++)
901    {
902      if((c=temp_array[small_i])!=0)
903      {
904        //res->idx_array[pos]=small_i;
905        //res->coef_array[pos]=temp_array[small_i];
906        (*(it_idx++))=small_i;
907        (*(it_coef++))=c;
908        //pos++;
909        non_zeros--;
910
911      }
912      if UNLIKELY(non_zeros==0) break;
913    }
914
915  }
916  ++it;
917}
918#endif
919return res;
920}
921template <class number_type> void add_coef_times_sparse(number_type* const temp_array,
922int temp_size,SparseRow<number_type>* row, number coef)
923{
924  int j;
925  number_type* const coef_array=row->coef_array;
926  int* const idx_array=row->idx_array;
927  const int len=row->len;
928  tgb_uint32 buffer[256];
929  const tgb_uint32 prime=n_GetChar(currRing->cf);
930  const tgb_uint32 c=F4mat_to_number_type(coef);
931  assume(!(n_IsZero(coef,currRing->cf)));
932  for(j=0;j<len;j=j+256)
933  {
934    const int bound=std::min(j+256,len);
935    int i;
936    int bpos=0;
937    for(i=j;i<bound;i++)
938    {
939      buffer[bpos++]=coef_array[i];
940    }
941    int bpos_bound=bound-j;
942    for(i=0;i<bpos_bound;i++)
943    {
944       buffer[i]*=c;
945     }
946    for(i=0;i<bpos_bound;i++)
947    {
948       buffer[i]=buffer[i]%prime;
949    }
950    bpos=0;
951    for(i=j;i<bound;i++)
952    {
953      int idx=idx_array[i];
954      assume(bpos<256);
955      assume(!(n_IsZero((number)(long) buffer[bpos],currRing->cf)));
956      temp_array[idx]=F4mat_to_number_type(n_Add((number)(long) temp_array[idx], (number)(long) buffer[bpos++],currRing->cf));
957      assume(idx<temp_size);
958    }
959
960  }
961}
962template <class number_type> void add_coef_times_dense(number_type* const temp_array,
963int temp_size,const number_type* row, int len,number coef)
964{
965  int j;
966  const number_type* const coef_array=row;
967  //int* const idx_array=row->idx_array;
968  //const int len=temp_size;
969  tgb_uint32 buffer[256];
970  const tgb_uint32 prime=n_GetChar(currRing->cf);
971  const tgb_uint32 c=F4mat_to_number_type(coef);
972  assume(!(n_IsZero(coef,currRing->cf)));
973  for(j=0;j<len;j=j+256)
974  {
975    const int bound=std::min(j+256,len);
976    int i;
977    int bpos=0;
978    for(i=j;i<bound;i++)
979    {
980      buffer[bpos++]=coef_array[i];
981    }
982    int bpos_bound=bound-j;
983    for(i=0;i<bpos_bound;i++)
984    {
985       buffer[i]*=c;
986     }
987    for(i=0;i<bpos_bound;i++)
988    {
989       buffer[i]=buffer[i]%prime;
990    }
991    bpos=0;
992    for(i=j;i<bound;i++)
993    {
994      //int idx=idx_array[i];
995      assume(bpos<256);
996      //assume(!(npIsZero((number) buffer[bpos])));
997      temp_array[i]=F4mat_to_number_type(n_Add((number)(long) temp_array[i], (number)(long) buffer[bpos++],currRing->cf));
998      assume(i<temp_size);
999    }
1000
1001  }
1002}
1003template <class number_type> void add_dense(number_type* const temp_array,
1004int temp_size,const number_type* row, int len)
1005{
1006  //int j;
1007  //const number_type* const coef_array=row;
1008  //int* const idx_array=row->idx_array;
1009  //const int len=temp_size;
1010  //tgb_uint32 buffer[256];
1011  //const tgb_uint32 prime=npPrimeM;
1012  //const tgb_uint32 c=F4mat_to_number_type(coef);
1013
1014  int i;
1015  for(i=0;i<len;i++)
1016  {
1017      temp_array[i]=F4mat_to_number_type(n_Add((number)(long) temp_array[i], (number)(long) row[i],currRing->cf));
1018      assume(i<temp_size);
1019  }
1020
1021}
1022template <class number_type> void sub_dense(number_type* const temp_array,
1023int temp_size,const number_type* row, int len)
1024{
1025  //int j;
1026  //const number_type* const coef_array=row;
1027  //int* const idx_array=row->idx_array;
1028  //const int len=temp_size;
1029  //tgb_uint32 buffer[256];
1030  //const tgb_uint32 prime=npPrimeM;
1031  //const tgb_uint32 c=F4mat_to_number_type(coef);
1032
1033  int i;
1034  for(i=0;i<len;i++)
1035  {
1036
1037      temp_array[i]=F4mat_to_number_type(n_Sub((number)(long) temp_array[i], (number)(long) row[i],currRing->cf));
1038      assume(i<temp_size);
1039  }
1040
1041}
1042
1043template <class number_type> void add_sparse(number_type* const temp_array,int temp_size,SparseRow<number_type>* row)
1044{
1045  int j;
1046
1047          number_type* const coef_array=row->coef_array;
1048          int* const idx_array=row->idx_array;
1049          const int len=row->len;
1050        for(j=0;j<len;j++)
1051        {
1052          int idx=idx_array[j];
1053          temp_array[idx]=F4mat_to_number_type(   (number_type)(long)n_Add((number) (long)temp_array[idx],(number)(long) coef_array[j],currRing->cf));
1054          assume(idx<temp_size);
1055        }
1056}
1057template <class number_type> void sub_sparse(number_type* const temp_array,int temp_size,SparseRow<number_type>* row)
1058{
1059  int j;
1060
1061          number_type* const coef_array=row->coef_array;
1062          int* const idx_array=row->idx_array;
1063          const int len=row->len;
1064        for(j=0;j<len;j++)
1065        {
1066          int idx=idx_array[j];
1067          temp_array[idx]=F4mat_to_number_type(  (number_type)(long) n_Sub((number) (long)temp_array[idx],(number)(long) coef_array[j],currRing->cf));
1068          assume(idx<temp_size);
1069        }
1070}
1071template <class number_type> SparseRow<number_type>* noro_red_to_non_poly_dense(MonRedResNP<number_type>* mon, int len,NoroCache<number_type>* cache)
1072{
1073  size_t temp_size_bytes=cache->nIrreducibleMonomials*sizeof(number_type)+8;//use 8bit int for testing
1074   assume(sizeof(int64)==8);
1075   cache->ensureTempBufferSize(temp_size_bytes);
1076   number_type* temp_array=(number_type*) cache->tempBuffer;//omalloc(cache->nIrreducibleMonomials*sizeof(number_type));
1077   int temp_size=cache->nIrreducibleMonomials;
1078   memset(temp_array,0,temp_size_bytes);
1079   number minus_one=n_Init(-1,currRing->cf);
1080   int i;
1081   for(i=0;i<len;i++)
1082   {
1083     MonRedResNP<number_type> red=mon[i];
1084     if ((red.ref))
1085     {
1086       if (red.ref->row)
1087       {
1088         SparseRow<number_type>* row=red.ref->row;
1089         number coef=red.coef;
1090         if (row->idx_array)
1091         {
1092           if (!((coef==(number)(long) 1)||(coef==minus_one)))
1093           {
1094             add_coef_times_sparse(temp_array,temp_size,row,coef);
1095           }
1096           else
1097           {
1098             if (coef==(number)(long) 1)
1099             {
1100               add_sparse(temp_array,temp_size,row);
1101             }
1102             else
1103             {
1104               sub_sparse(temp_array,temp_size,row);
1105             }
1106           }
1107         }
1108         else
1109         //TODO: treat, 1,-1
1110         if (!((coef==(number)(long) 1)||(coef==minus_one)))
1111         {
1112          add_coef_times_dense(temp_array,temp_size,row->coef_array,row->len,coef);
1113         }
1114         else
1115         {
1116           if (coef==(number)(long)1)
1117             add_dense(temp_array,temp_size,row->coef_array,row->len);
1118           else
1119           {
1120             assume(coef==minus_one);
1121             sub_dense(temp_array,temp_size,row->coef_array,row->len);
1122             //add_coef_times_dense(temp_array,temp_size,row->coef_array,row->len,coef);
1123           }
1124         }
1125       }
1126       else
1127       {
1128         if (red.ref->value_len==NoroCache<number_type>::backLinkCode)
1129         {
1130           temp_array[red.ref->term_index]=F4mat_to_number_type( n_Add((number)(long) temp_array[red.ref->term_index],red.coef,currRing->cf));
1131         }
1132         else
1133         {
1134           //PrintS("third case\n");
1135         }
1136       }
1137     }
1138   }
1139   int non_zeros=0;
1140   for(i=0;i<cache->nIrreducibleMonomials;i++)
1141   {
1142     //if (!(temp_array[i]==0)){
1143     //  non_zeros++;
1144     //}
1145     assume(((temp_array[i]!=0)==0)|| (((temp_array[i]!=0)==1)));
1146     non_zeros+=(temp_array[i]!=0);
1147   }
1148
1149   if (non_zeros==0)
1150   {
1151     //omfree(mon);
1152     return NULL;
1153   }
1154   SparseRow<number_type>* res=new SparseRow<number_type>(temp_size,temp_array);//convert_to_sparse_row(temp_array,temp_size, non_zeros);
1155
1156   //omfree(temp_array);
1157
1158
1159   return res;
1160}
1161template<class number_type> class CoefIdx
1162{
1163public:
1164  number_type coef;
1165  int idx;
1166  bool operator<(const CoefIdx<number_type>& other) const
1167  {
1168    return (idx<other.idx);
1169  }
1170};
1171template<class number_type> void write_coef_times_xx_idx_to_buffer(CoefIdx<number_type>* const pairs,int& pos,int* const idx_array, number_type* const coef_array,const int rlen, const number coef)
1172{
1173  int j;
1174  for(j=0;j<rlen;j++)
1175  {
1176    assume(coef_array[j]!=0);
1177    CoefIdx<number_type> ci;
1178    ci.coef=F4mat_to_number_type(n_Mult((number)(long) coef,(number)(long) coef_array[j],currRing->cf));
1179    ci.idx=idx_array[j];
1180    pairs[pos++]=ci;
1181  }
1182}
1183template<class number_type> void write_coef_times_xx_idx_to_buffer_dense(CoefIdx<number_type>* const pairs,int& pos, number_type* const coef_array,const int rlen, const number coef)
1184{
1185  int j;
1186
1187  for(j=0;j<rlen;j++)
1188  {
1189    if (coef_array[j]!=0)
1190    {
1191      assume(coef_array[j]!=0);
1192      CoefIdx<number_type> ci;
1193      ci.coef=F4mat_to_number_type(n_Mult((number)(long) coef,(number)(long) coef_array[j],currRing->cf));
1194      assume(ci.coef!=0);
1195      ci.idx=j;
1196      pairs[pos++]=ci;
1197    }
1198  }
1199}
1200template<class number_type> void write_coef_idx_to_buffer_dense(CoefIdx<number_type>* const pairs,int& pos, number_type* const coef_array,const int rlen)
1201{
1202  int j;
1203
1204  for(j=0;j<rlen;j++)
1205  {
1206    if (coef_array[j]!=0)
1207    {
1208      assume(coef_array[j]!=0);
1209      CoefIdx<number_type> ci;
1210      ci.coef=coef_array[j];
1211      assume(ci.coef!=0);
1212      ci.idx=j;
1213      pairs[pos++]=ci;
1214    }
1215  }
1216}
1217
1218template<class number_type> void write_minus_coef_idx_to_buffer_dense(CoefIdx<number_type>* const pairs,int& pos, number_type* const coef_array,const int rlen)
1219{
1220  int j;
1221
1222  for(j=0;j<rlen;j++)
1223  {
1224    if (coef_array[j]!=0)
1225    {
1226      assume(coef_array[j]!=0);
1227      CoefIdx<number_type> ci;
1228      ci.coef=F4mat_to_number_type(n_Neg((number)(long) coef_array[j],currRing->cf));
1229      assume(ci.coef!=0);
1230      ci.idx=j;
1231      pairs[pos++]=ci;
1232    }
1233  }
1234}
1235template<class number_type> void write_coef_idx_to_buffer(CoefIdx<number_type>* const pairs,int& pos,int* const idx_array, number_type* const coef_array,const int rlen)
1236{
1237  int j;
1238  for(j=0;j<rlen;j++)
1239  {
1240    assume(coef_array[j]!=0);
1241    CoefIdx<number_type> ci;
1242    ci.coef=coef_array[j];
1243    ci.idx=idx_array[j];
1244    pairs[pos++]=ci;
1245  }
1246}
1247
1248template<class number_type> void write_minus_coef_idx_to_buffer(CoefIdx<number_type>* const pairs,int& pos,int* const idx_array, number_type* const coef_array,const int rlen)
1249{
1250  int j;
1251  for(j=0;j<rlen;j++)
1252  {
1253    assume(coef_array[j]!=0);
1254    CoefIdx<number_type> ci;
1255    ci.coef=F4mat_to_number_type(n_Neg((number)coef_array[j],currRing->cf));
1256    ci.idx=idx_array[j];
1257    pairs[pos++]=ci;
1258  }
1259}
1260template <class number_type> SparseRow<number_type>* noro_red_to_non_poly_sparse(MonRedResNP<number_type>* mon, int len,NoroCache<number_type>* cache)
1261{
1262  int i;
1263  int together=0;
1264  for(i=0;i<len;i++)
1265  {
1266    MonRedResNP<number_type> red=mon[i];
1267    if ((red.ref) &&( red.ref->row))
1268    {
1269      together+=red.ref->row->len;
1270    }
1271    else
1272    {
1273      if ((red.ref) &&(red.ref->value_len==NoroCache<number_type>::backLinkCode))
1274      together++;
1275    }
1276  }
1277  //PrintS("here\n");
1278  if (together==0) return 0;
1279  //PrintS("there\n");
1280  cache->ensureTempBufferSize(together*sizeof(CoefIdx<number_type>));
1281  CoefIdx<number_type>* pairs=(CoefIdx<number_type>*) cache->tempBuffer; //omalloc(together*sizeof(CoefIdx<number_type>));
1282  int pos=0;
1283  const number one=n_Init(1, currRing->cf);
1284  const number minus_one=n_Init(-1, currRing->cf);
1285  for(i=0;i<len;i++)
1286  {
1287    MonRedResNP<number_type> red=mon[i];
1288    if ((red.ref) &&( red.ref->row))
1289    {
1290      //together+=red.ref->row->len;
1291      int* idx_array=red.ref->row->idx_array;
1292      number_type* coef_array=red.ref->row->coef_array;
1293      int rlen=red.ref->row->len;
1294      number coef=red.coef;
1295      if (idx_array)
1296      {
1297        if ((coef!=one)&&(coef!=minus_one))
1298        {
1299          write_coef_times_xx_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen, coef);
1300        }
1301        else
1302        {
1303          if (coef==one)
1304          {
1305            write_coef_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen);
1306          }
1307          else
1308          {
1309            assume(coef==minus_one);
1310            write_minus_coef_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen);
1311          }
1312        }
1313      }
1314      else
1315      {
1316        if ((coef!=one)&&(coef!=minus_one))
1317        {
1318          write_coef_times_xx_idx_to_buffer_dense(pairs,pos,coef_array,rlen,coef);
1319        }
1320        else
1321        {
1322          if (coef==one)
1323            write_coef_idx_to_buffer_dense(pairs,pos,coef_array,rlen);
1324          else
1325          {
1326            assume(coef==minus_one);
1327            write_minus_coef_idx_to_buffer_dense(pairs,pos,coef_array,rlen);
1328          }
1329        }
1330      }
1331    }
1332    else
1333    {
1334      if ((red.ref) &&(red.ref->value_len==NoroCache<number_type>::backLinkCode))
1335      {
1336        CoefIdx<number_type> ci;
1337        ci.coef=F4mat_to_number_type(red.coef);
1338        ci.idx=red.ref->term_index;
1339        pairs[pos++]=ci;
1340      }
1341    }
1342  }
1343  assume(pos<=together);
1344  together=pos;
1345
1346  std::sort(pairs,pairs+together);
1347
1348  int act=0;
1349
1350  assume(pairs[0].coef!=0);
1351  for(i=1;i<together;i++)
1352  {
1353    if (pairs[i].idx!=pairs[act].idx)
1354    {
1355      if (pairs[act].coef!=0)
1356      {
1357        act=act+1;
1358      }
1359      pairs[act]=pairs[i];
1360    }
1361    else
1362    {
1363      pairs[act].coef=F4mat_to_number_type(n_Add((number)(long)pairs[act].coef,(number)(long)pairs[i].coef,currRing->cf));
1364    }
1365  }
1366
1367  if (pairs[act].coef==0)
1368  {
1369    act--;
1370  }
1371  int sparse_row_len=act+1;
1372  //Print("res len:%d",sparse_row_len);
1373  if (sparse_row_len==0) {return NULL;}
1374  SparseRow<number_type>* res=new SparseRow<number_type>(sparse_row_len);
1375  {
1376    number_type* coef_array=res->coef_array;
1377    int* idx_array=res->idx_array;
1378    for(i=0;i<sparse_row_len;i++)
1379    {
1380      idx_array[i]=pairs[i].idx;
1381      coef_array[i]=pairs[i].coef;
1382    }
1383  }
1384  //omfree(pairs);
1385
1386  return res;
1387}
1388template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c){
1389  assume(len==pLength(p));
1390  if (p==NULL)
1391  {
1392    len=0;
1393    return NULL;
1394  }
1395
1396  MonRedResNP<number_type>* mon=(MonRedResNP<number_type>*) omalloc(len*sizeof(MonRedResNP<number_type>));
1397  int i=0;
1398  double max_density=0.0;
1399  while(p!=NULL)
1400  {
1401    poly t=p;
1402    pIter(p);
1403    pNext(t)=NULL;
1404
1405    MonRedResNP<number_type> red=noro_red_mon_to_non_poly(t,cache,c);
1406    if ((red.ref) && (red.ref->row))
1407    {
1408      double act_density=(double) red.ref->row->len;
1409      act_density/=(double) cache->nIrreducibleMonomials;
1410      max_density=std::max(act_density,max_density);
1411    }
1412    mon[i]=red;
1413    i++;
1414  }
1415
1416  assume(i==len);
1417  len=i;
1418  bool dense=true;
1419  if (max_density<0.3) dense=false;
1420  if (dense){
1421    SparseRow<number_type>* res=noro_red_to_non_poly_dense(mon,len,cache);
1422    omfree(mon);
1423    return res;
1424  } else   {
1425      SparseRow<number_type>* res=noro_red_to_non_poly_sparse(mon,len,cache);
1426      omfree(mon);
1427      return res;
1428    }
1429  //in the loop before nIrreducibleMonomials increases, so position here is important
1430
1431}
1432#endif
1433static wlen_type pair_weighted_length(int i, int j, slimgb_alg* c);
1434wlen_type pELength(poly p, ring r);
1435int terms_sort_crit(const void* a, const void* b);
1436//void simplest_gauss_modp(number* a, int nrows,int ncols);
1437// a: a[0,0],a[0,1]....a[nrows-1,ncols-1]
1438// assume: field is Zp
1439#ifdef USE_NORO
1440
1441
1442template <class number_type > void write_poly_to_row(number_type* row, poly h, poly*terms, int tn, ring r){
1443  //poly* base=row;
1444  while(h!=NULL){
1445    //Print("h:%i\n",h);
1446    number coef=p_GetCoeff(h,r);
1447    poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit);
1448    assume(ptr_to_h!=NULL);
1449    int pos=ptr_to_h-terms;
1450    row[pos]=F4mat_to_number_type(coef);
1451    //number_type_array[base+pos]=coef;
1452    pIter(h);
1453  }
1454}
1455template <class number_type > poly row_to_poly(number_type* row, poly* terms, int tn, ring r){
1456  poly h=NULL;
1457  int j;
1458  number_type zero=0;//;npInit(0);
1459  for(j=tn-1;j>=0;j--){
1460    if (!(zero==(row[j]))){
1461      poly t=terms[j];
1462      t=p_LmInit(t,r);
1463      p_SetCoeff(t,(number)(long) row[j],r);
1464      pNext(t)=h;
1465      h=t;
1466    }
1467
1468  }
1469  return h;
1470}
1471template <class number_type > int modP_lastIndexRow(number_type* row,int ncols)
1472{
1473  int lastIndex;
1474  const number_type zero=0;//npInit(0);
1475  for(lastIndex=ncols-1;lastIndex>=0;lastIndex--)
1476  {
1477    if (!(row[lastIndex]==zero))
1478    {
1479      return lastIndex;
1480    }
1481  }
1482  return -1;
1483}
1484template <class number_type> int term_nodes_sort_crit(const void* a, const void* b)
1485{
1486  return -pLmCmp(((TermNoroDataNode<number_type>*) a)->t,((TermNoroDataNode<number_type>*) b)->t);
1487}
1488
1489template <class number_type>class ModPMatrixBackSubstProxyOnArray;
1490template <class number_type > class ModPMatrixProxyOnArray
1491{
1492public:
1493  friend class ModPMatrixBackSubstProxyOnArray<number_type>;
1494
1495  int ncols,nrows;
1496  ModPMatrixProxyOnArray(number_type* array, int nnrows, int nncols)
1497  {
1498    this->ncols=nncols;
1499    this->nrows=nnrows;
1500    rows=(number_type**) omalloc(nnrows*sizeof(number_type*));
1501    startIndices=(int*)omalloc(nnrows*sizeof(int));
1502    int i;
1503    for(i=0;i<nnrows;i++)
1504    {
1505      rows[i]=array+(i*nncols);
1506      updateStartIndex(i,-1);
1507    }
1508  }
1509  ~ModPMatrixProxyOnArray()
1510  {
1511    omfree(rows);
1512    omfree(startIndices);
1513  }
1514
1515  void permRows(int i, int j)
1516  {
1517    number_type* h=rows[i];
1518    rows[i]=rows[j];
1519    rows[j]=h;
1520    int hs=startIndices[i];
1521    startIndices[i]=startIndices[j];
1522    startIndices[j]=hs;
1523  }
1524  void multiplyRow(int row, number_type coef)
1525  {
1526    int i;
1527    number_type* row_array=rows[row];
1528    for(i=startIndices[row];i<ncols;i++)
1529    {
1530      row_array[i]=F4mat_to_number_type(n_Mult((number)(long) row_array[i],(number)(long) coef,currRing->cf));
1531    }
1532  }
1533  void reduceOtherRowsForward(int r)
1534  {
1535    //assume rows "under r" have bigger or equal start index
1536    number_type* row_array=rows[r];
1537    number_type zero=F4mat_to_number_type((number)0 /*npInit(0, currRing)*/);
1538    int start=startIndices[r];
1539    number_type coef=row_array[start];
1540    assume(start<ncols);
1541    int other_row;
1542    assume(!(n_IsZero((number)(long) row_array[start],currRing->cf)));
1543    if (!(n_IsOne((number)(long) coef,currRing->cf)))
1544      multiplyRow(r,F4mat_to_number_type(n_Invers((number)(long) coef,currRing->cf)));
1545    assume(n_IsOne((number)(long) row_array[start],currRing->cf));
1546    int lastIndex=modP_lastIndexRow(row_array, ncols);
1547    number minus_one=n_Init(-1, currRing->cf);
1548    for (other_row=r+1;other_row<nrows;other_row++)
1549    {
1550      assume(startIndices[other_row]>=start);
1551      if (startIndices[other_row]==start)
1552      {
1553        int i;
1554        number_type* other_row_array=rows[other_row];
1555        number coef2=n_Neg((number)(long) other_row_array[start],currRing->cf);
1556        if (coef2==minus_one)
1557        {
1558          for(i=start;i<=lastIndex;i++)
1559          {
1560            if (row_array[i]!=zero)
1561              other_row_array[i]=F4mat_to_number_type(n_Sub((number)(long) other_row_array[i], (number)(long) row_array[i],currRing->cf));
1562          }
1563      }
1564      else
1565      {
1566          //assume(FALSE);
1567          for(i=start;i<=lastIndex;i++)
1568          {
1569            if (row_array[i]!=zero)
1570            other_row_array[i]=F4mat_to_number_type(n_Add(n_Mult(coef2,(number)(long) row_array[i],currRing->cf),(number)(long) other_row_array[i],currRing->cf));
1571          }
1572        }
1573        updateStartIndex(other_row,start);
1574        assume(n_IsZero((number)(long) other_row_array[start],currRing->cf));
1575      }
1576    }
1577  }
1578  void updateStartIndex(int row,int lower_bound)
1579  {
1580    number_type* row_array=rows[row];
1581    assume((lower_bound<0)||(n_IsZero((number)(long) row_array[lower_bound],currRing->cf)));
1582    int i;
1583    //number_type zero=npInit(0);
1584    for(i=lower_bound+1;i<ncols;i++)
1585    {
1586      if (!(row_array[i]==0))
1587        break;
1588    }
1589    startIndices[row]=i;
1590  }
1591  int getStartIndex(int row)
1592  {
1593    return startIndices[row];
1594  }
1595  BOOLEAN findPivot(int &r, int &c)
1596  {
1597    //row>=r, col>=c
1598
1599    while(c<ncols)
1600    {
1601      int i;
1602      for(i=r;i<nrows;i++)
1603      {
1604        assume(startIndices[i]>=c);
1605        if (startIndices[i]==c)
1606        {
1607          //r=i;
1608          if (r!=i)
1609            permRows(r,i);
1610          return TRUE;
1611        }
1612      }
1613      c++;
1614    }
1615    return FALSE;
1616  }
1617protected:
1618  number_type** rows;
1619  int* startIndices;
1620};
1621template <class number_type > class ModPMatrixBackSubstProxyOnArray
1622{
1623  int *startIndices;
1624  number_type** rows;
1625  int *lastReducibleIndices;
1626  int ncols;
1627  int nrows;
1628  int nonZeroUntil;
1629public:
1630  void multiplyRow(int row, number_type coef)
1631  {
1632    int i;
1633    number_type* row_array=rows[row];
1634    for(i=startIndices[row];i<ncols;i++)
1635    {
1636      row_array[i]=F4mat_to_number_type(n_Mult((number)(long) row_array[i],(number)(long) coef,currRing->cf));
1637    }
1638  }
1639  ModPMatrixBackSubstProxyOnArray<number_type> (ModPMatrixProxyOnArray<number_type> & p)
1640  {
1641//  (number_type* array, int nrows, int ncols, int* startIndices, number_type** rows){
1642    //we borrow some parameters ;-)
1643    //we assume, that nobody changes the order of the rows
1644    this->startIndices=p.startIndices;
1645    this->rows=p.rows;
1646    this->ncols=p.ncols;
1647    this->nrows=p.nrows;
1648    lastReducibleIndices=(int*) omalloc(nrows*sizeof(int));
1649    nonZeroUntil=0;
1650    while(nonZeroUntil<nrows)
1651    {
1652      if (startIndices[nonZeroUntil]<ncols)
1653      {
1654        nonZeroUntil++;
1655      }
1656      else break;
1657    }
1658    if (TEST_OPT_PROT)
1659      Print("rank:%i\n",nonZeroUntil);
1660    nonZeroUntil--;
1661    int i;
1662    for(i=0;i<=nonZeroUntil;i++)
1663    {
1664      assume(startIndices[i]<ncols);
1665      assume(!(n_IsZero((number)(long) rows[i][startIndices[i]],currRing->cf)));
1666      assume(startIndices[i]>=i);
1667      updateLastReducibleIndex(i,nonZeroUntil+1);
1668    }
1669  }
1670  void updateLastReducibleIndex(int r, int upper_bound)
1671  {
1672    number_type* row_array=rows[r];
1673    if (upper_bound>nonZeroUntil) upper_bound=nonZeroUntil+1;
1674    int i;
1675    const number_type zero=0;//npInit(0);
1676    for(i=upper_bound-1;i>r;i--)
1677    {
1678      int start=startIndices[i];
1679      assume(start<ncols);
1680      if (!(row_array[start]==zero))
1681      {
1682        lastReducibleIndices[r]=start;
1683        return;
1684      }
1685    }
1686    lastReducibleIndices[r]=-1;
1687  }
1688  void backwardSubstitute(int r)
1689  {
1690    int start=startIndices[r];
1691    assume(start<ncols);
1692    number_type zero=0;//npInit(0);
1693    number_type* row_array=rows[r];
1694    assume((!(n_IsZero((number)(long) row_array[start],currRing->cf))));
1695    assume(start<ncols);
1696    int other_row;
1697    if (!(n_IsOne((number)(long) row_array[r],currRing->cf)))
1698    {
1699      //it should be one, but this safety is not expensive
1700      multiplyRow(r, F4mat_to_number_type(n_Invers((number)(long) row_array[start],currRing->cf)));
1701    }
1702    int lastIndex=modP_lastIndexRow(row_array, ncols);
1703    assume(lastIndex<ncols);
1704    assume(lastIndex>=0);
1705    for(other_row=r-1;other_row>=0;other_row--)
1706    {
1707      assume(lastReducibleIndices[other_row]<=start);
1708      if (lastReducibleIndices[other_row]==start)
1709      {
1710        number_type* other_row_array=rows[other_row];
1711        number coef=n_Neg((number)(long) other_row_array[start],currRing->cf);
1712        assume(!(n_IsZero(coef,currRing->cf)));
1713        int i;
1714        assume(start>startIndices[other_row]);
1715        for(i=start;i<=lastIndex;i++)
1716        {
1717          if (row_array[i]!=zero)
1718            other_row_array[i]=F4mat_to_number_type(n_Add(n_Mult(coef,(number)(long)row_array[i],currRing->cf),(number)(long)other_row_array[i],currRing->cf));
1719        }
1720        updateLastReducibleIndex(other_row,r);
1721      }
1722    }
1723  }
1724  ~ModPMatrixBackSubstProxyOnArray<number_type>()
1725  {
1726    omfree(lastReducibleIndices);
1727  }
1728  void backwardSubstitute()
1729  {
1730    int i;
1731    for(i=nonZeroUntil;i>0;i--)
1732    {
1733      backwardSubstitute(i);
1734    }
1735  }
1736};
1737template <class number_type > void simplest_gauss_modp(number_type* a, int nrows,int ncols)
1738{
1739  //use memmoves for changing rows
1740  //if (TEST_OPT_PROT)
1741  //    PrintS("StartGauss\n");
1742  ModPMatrixProxyOnArray<number_type> mat(a,nrows,ncols);
1743
1744  int c=0;
1745  int r=0;
1746  while(mat.findPivot(r,c)){
1747    //int pivot=find_pivot()
1748      mat.reduceOtherRowsForward(r);
1749    r++;
1750    c++;
1751  }
1752  ModPMatrixBackSubstProxyOnArray<number_type> backmat(mat);
1753  backmat.backwardSubstitute();
1754  //backward substitutions
1755  //if (TEST_OPT_PROT)
1756  //PrintS("StopGauss\n");
1757}
1758//int term_nodes_sort_crit(const void* a, const void* b);
1759template <class number_type> void noro_step(poly*p,int &pn,slimgb_alg* c){
1760  //Print("Input rows %d\n",pn);
1761  int j;
1762  if (TEST_OPT_PROT)
1763  {
1764    Print("Input rows %d\n",pn);
1765  }
1766
1767  NoroCache<number_type> cache;
1768
1769  SparseRow<number_type> ** srows=(SparseRow<number_type>**) omAlloc(pn*sizeof(SparseRow<number_type>*));
1770  int non_zeros=0;
1771  for(j=0;j<pn;j++)
1772  {
1773    poly h=p[j];
1774    int h_len=pLength(h);
1775    //number coef;
1776    srows[non_zeros]=noro_red_to_non_poly_t<number_type>(h,h_len,&cache,c);
1777    if (srows[non_zeros]!=NULL) non_zeros++;
1778  }
1779  std::vector<DataNoroCacheNode<number_type>*> irr_nodes;
1780  cache.collectIrreducibleMonomials(irr_nodes);
1781  //now can build up terms array
1782  //Print("historic irred Mon%d\n",cache.nIrreducibleMonomials);
1783  int n=irr_nodes.size();//cache.countIrreducibleMonomials();
1784  cache.nIrreducibleMonomials=n;
1785  if (TEST_OPT_PROT)
1786  {
1787    Print("Irred Mon:%d\n",n);
1788    Print("red Mon:%d\n",cache.nReducibleMonomials);
1789  }
1790  TermNoroDataNode<number_type>* term_nodes=(TermNoroDataNode<number_type>*) omalloc(n*sizeof(TermNoroDataNode<number_type>));
1791
1792  for(j=0;j<n;j++)
1793  {
1794    assume(irr_nodes[j]!=NULL);
1795    assume(irr_nodes[j]->value_len==NoroCache<number_type>::backLinkCode);
1796    term_nodes[j].t=irr_nodes[j]->value_poly;
1797    assume(term_nodes[j].t!=NULL);
1798    term_nodes[j].node=irr_nodes[j];
1799  }
1800
1801  qsort(term_nodes,n,sizeof(TermNoroDataNode<number_type>),term_nodes_sort_crit<number_type>);
1802  poly* terms=(poly*) omalloc(n*sizeof(poly));
1803
1804  int* old_to_new_indices=(int*) omalloc(cache.nIrreducibleMonomials*sizeof(int));
1805  for(j=0;j<n;j++)
1806  {
1807    old_to_new_indices[term_nodes[j].node->term_index]=j;
1808    term_nodes[j].node->term_index=j;
1809    terms[j]=term_nodes[j].t;
1810  }
1811
1812  //if (TEST_OPT_PROT)
1813  //  Print("Evaluate Rows \n");
1814  pn=non_zeros;
1815  number_type* number_array=(number_type*) omalloc0(n*pn*sizeof(number_type));
1816
1817  for(j=0;j<pn;j++)
1818  {
1819    int i;
1820    number_type* row=number_array+n*j;
1821    /*for(i=0;i<n;i++)
1822    {
1823      row[i]=zero;
1824    }*/
1825
1826    SparseRow<number_type>* srow=srows[j];
1827
1828    if (srow)
1829    {
1830      int* const idx_array=srow->idx_array;
1831      number_type* const coef_array=srow->coef_array;
1832      const int len=srow->len;
1833      if (srow->idx_array)
1834      {
1835        for(i=0;i<len;i++)
1836        {
1837         int idx=old_to_new_indices[idx_array[i]];
1838         row[idx]=F4mat_to_number_type(coef_array[i]);
1839        }
1840      }
1841      else
1842      {
1843        for(i=0;i<len;i++)
1844        {
1845          row[old_to_new_indices[i]]=F4mat_to_number_type(coef_array[i]);
1846        }
1847      }
1848      delete srow;
1849    }
1850  }
1851
1852  //static int export_n=0;
1853  //export_mat(number_array,pn,n,"mat%i.py",++export_n);
1854  simplest_gauss_modp(number_array,pn,n);
1855
1856  int p_pos=0;
1857  for(j=0;j<pn;j++){
1858    poly h=row_to_poly(number_array+j*n,terms,n,c->r);
1859    if(h!=NULL){
1860      p[p_pos++]=h;
1861    }
1862  }
1863  pn=p_pos;
1864  omfree(terms);
1865  omfree(term_nodes);
1866  omfree(number_array);
1867  #ifdef NORO_NON_POLY
1868  omfree(srows);
1869  omfree(old_to_new_indices);
1870  #endif
1871  //don't forget the rank
1872
1873}
1874
1875template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type> *>& res){
1876  int i;
1877  for(i=0;i<root.branches_len;i++){
1878    collectIrreducibleMonomials(1,root.branches[i],res);
1879  }
1880}
1881template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials(int level, NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>*>& res){
1882  assume(level>=0);
1883  if (node==NULL) return;
1884  if (level<(currRing->N))
1885  {
1886    int i;
1887    for(i=0;i<node->branches_len;i++)
1888    {
1889      collectIrreducibleMonomials(level+1,node->branches[i],res);
1890    }
1891  }
1892  else
1893  {
1894    DataNoroCacheNode<number_type>* dn=(DataNoroCacheNode<number_type>*) node;
1895    if (dn->value_len==backLinkCode)
1896    {
1897      res.push_back(dn);
1898    }
1899  }
1900}
1901
1902template<class number_type> DataNoroCacheNode<number_type>* NoroCache<number_type>::getCacheReference(poly term){
1903  int i;
1904  NoroCacheNode* parent=&root;
1905  for(i=1;i<(currRing->N);i++){
1906    parent=parent->getBranch(p_GetExp(term,i,currRing));
1907    if (!(parent)){
1908      return NULL;
1909    }
1910  }
1911  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1912  return res_holder;
1913}
1914template<class number_type> poly NoroCache<number_type>::lookup(poly term, BOOLEAN& succ, int & len){
1915  int i;
1916  NoroCacheNode* parent=&root;
1917  for(i=1;i<(currRing->N);i++){
1918    parent=parent->getBranch(p_GetExp(term,i,currRing));
1919    if (!(parent)){
1920      succ=FALSE;
1921      return NULL;
1922    }
1923  }
1924  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1925  if (res_holder){
1926    succ=TRUE;
1927    if ((res_holder->value_len==backLinkCode)){
1928      len=1;
1929      return term;
1930    }
1931    len=res_holder->value_len;
1932    return res_holder->value_poly;
1933  } else {
1934    succ=FALSE;
1935    return NULL;
1936  }
1937}
1938#endif
1939
1940#endif
Note: See TracBrowser for help on using the repository browser.