source: git/kernel/tgb_internal.h @ a5b80a

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