source: git/kernel/tgb_internal.h @ 1b2d84

jengelh-datetimespielwiese
Last change on this file since 1b2d84 was 1b2d84, checked in by Michael Brickenstein <bricken@…>, 16 years ago
*bricken: si_max ->std::max git-svn-id: file:///usr/local/Singular/svn/trunk@9903 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 42.0 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: tgb_internal.h,v 1.66 2007-02-26 15:14:47 bricken Exp $ */
8/*
9 * ABSTRACT: tgb internal .h file
10*/
11#include <omalloc.h>
12#include "p_polys.h"
13
14#include "ideals.h"
15#include "ring.h"
16#include "febase.h"
17#include "structs.h"
18#include "polys.h"
19#include "stdlib.h"
20#include <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 "boost/dynamic_bitset.hpp"
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 "kutil.h"
50#include "kInline.cc"
51#include "kstd1.h"
52#include "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    impl=p;
76  }
77  PolySimple(){
78    impl=NULL;
79  }
80  PolySimple(const PolySimple& a){
81    //impl=p_Copy(a.impl,currRing);
82    impl=a.impl;
83  }
84  PolySimple& operator=(const PolySimple& p2){
85    //p_Delete(&impl,currRing);
86    //impl=p_Copy(p2.impl,currRing);
87    impl=p2.impl;
88    return *this;
89  }
90  ~PolySimple(){
91    //p_Delete(&impl,currRing);
92  }
93  bool operator< (const PolySimple& other) const{
94    return pLmCmp(impl,other.impl)<0;
95  }
96  bool operator==(const PolySimple& other){
97    return pLmEqual(impl,other.impl);
98  }
99  poly impl;
100 
101};
102template<class number_type> class DataNoroCacheNode;
103/*class MonRedRes{
104public:
105  poly p;
106  number coef;
107  BOOLEAN changed;
108  int len;
109  BOOLEAN onlyBorrowed;
110  bool operator<(const MonRedRes& other) const{
111    int cmp=p_LmCmp(p,other.p,currRing);
112    if ((cmp<0)||((cmp==0)&&((onlyBorrowed)&&(!(other.onlyBorrowed))))){
113      return true;
114    } else return false;
115  }
116  DataNoroCacheNode* ref;
117  MonRedRes(){
118    ref=NULL;
119    p=NULL;
120  }
121};*/
122template <class number_type> class MonRedResNP{
123public:
124  number coef;
125
126
127  DataNoroCacheNode<number_type>* ref;
128  MonRedResNP(){
129    ref=NULL;
130  }
131};
132struct sorted_pair_node{
133  //criterium, which is stable 0. small lcm 1. small i 2. small j
134  wlen_type expected_length;
135  poly lcm_of_lm;
136  int i;
137  int j;
138  int deg;
139 
140 
141};
142#ifdef NORO_CACHE
143#ifndef NORO_NON_POLY
144class NoroPlaceHolder{
145public:
146  DataNoroCacheNode* ref;
147  number coef;
148};
149#endif
150#endif
151//static ideal debug_Ideal;
152
153
154struct poly_list_node{
155  poly p;
156  poly_list_node* next;
157};
158
159struct int_pair_node{
160  int_pair_node* next;
161  int a;
162  int b;
163};
164struct monom_poly{
165  poly m;
166  poly f;
167};
168struct mp_array_list{
169  monom_poly* mp;
170  int size;
171  mp_array_list* next;
172};
173
174
175struct poly_array_list{
176  poly* p;
177  int size;
178  poly_array_list* next;
179};
180class slimgb_alg
181{
182  public:
183    slimgb_alg(ideal I, int syz_comp,BOOLEAN F4);
184                void introduceDelayedPairs(poly* pa,int s);
185    virtual ~slimgb_alg();
186    void cleanDegs(int lower, int upper);
187#ifndef HAVE_BOOST
188#ifdef USE_STDVECBOOL
189  vector<vector<bool> > states;
190#else
191  char** states;
192#endif
193#else
194  vector<dynamic_bitset<> > states;
195#endif
196  ideal add_later;
197  ideal S;
198  ring r;
199  int* lengths;
200  wlen_type* weighted_lengths;
201  long* short_Exps;
202  kStrategy strat;
203  int* T_deg;
204  int* T_deg_full;
205  poly tmp_lm;
206  poly* tmp_pair_lm;
207  sorted_pair_node** tmp_spn;
208  poly* expandS;
209  poly* gcd_of_terms;
210  int_pair_node* soon_free;
211  sorted_pair_node** apairs;
212  #if 0
213  BOOLEAN* modifiedS;
214  #endif
215  #ifdef TGB_RESORT_PAIRS
216  bool* replaced;
217  #endif
218  poly_list_node* to_destroy;
219  //for F4
220  mp_array_list* F;
221  poly_array_list* F_minus;
222
223  //end for F4
224#ifdef HEAD_BIN
225  struct omBin_s*   HeadBin;
226#endif
227  unsigned int reduction_steps;
228  int n;
229  //! array_lengths should be greater equal n;
230  int syz_comp;
231  int array_lengths; 
232  int normal_forms;
233  int current_degree;
234  int Rcounter;
235  int last_index;
236  int max_pairs;
237  int pair_top;
238  int easy_product_crit;
239  int extended_product_crit;
240  int average_length;
241  int lastDpBlockStart;
242  int lastCleanedDeg;
243  BOOLEAN isDifficultField;
244  BOOLEAN completed;
245  BOOLEAN is_homog;
246  BOOLEAN tailReductions;
247  BOOLEAN eliminationProblem;
248  BOOLEAN F4_mode;
249  BOOLEAN nc;
250  #ifdef TGB_RESORT_PAIRS
251  BOOLEAN used_b;
252  #endif
253};
254class red_object{
255 public:
256  kBucket_pt bucket;
257  poly p;
258  unsigned long sev;
259  int sugar;
260  void flatten();
261  void validate();
262  wlen_type initial_quality;
263  void adjust_coefs(number c_r, number c_ac_r);
264  wlen_type guess_quality(slimgb_alg* c);
265  int clear_to_poly();
266  void canonicalize();
267};
268
269
270enum calc_state
271  {
272    UNCALCULATED,
273    HASTREP//,
274    //UNIMPORTANT,
275    //SOONTREP
276  };
277static BOOLEAN pair_cmp(sorted_pair_node* a,sorted_pair_node* b);
278template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set);
279static int add_to_reductors(slimgb_alg* c, poly h, int len, int ecart, BOOLEAN simplified=FALSE);
280static int bucket_guess(kBucket* bucket);
281static poly redNFTail (poly h,const int sl,kStrategy strat, int len);
282static poly redNF2 (poly h,slimgb_alg* c , int &len, number&  m,int n=0);
283void free_sorted_pair_node(sorted_pair_node* s, ring r);
284static void shorten_tails(slimgb_alg* c, poly monom);
285static void replace_pair(int & i, int & j, slimgb_alg* c);
286//static sorted_pair_node** add_to_basis(poly h, int i, int j,slimgb_alg* c, int* ip=NULL);
287static void do_this_spoly_stuff(int i,int j,slimgb_alg* c);
288//ideal t_rep_gb(ring r,ideal arg_I);
289static BOOLEAN has_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* state);
290static int* make_connections(int from, poly bound, slimgb_alg* c);
291static int* make_connections(int from, int to, poly bound, slimgb_alg* c);
292void now_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
293static void soon_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
294static int pLcmDeg(poly a, poly b);
295static int simple_posInS (kStrategy strat, poly p,int len, wlen_type wlen);
296static BOOLEAN find_next_pair(slimgb_alg* c, BOOLEAN go_higher=TRUE);
297
298static sorted_pair_node* pop_pair(slimgb_alg* c);
299static BOOLEAN no_pairs(slimgb_alg* c);
300void clean_top_of_pair_list(slimgb_alg* c);
301static void super_clean_top_of_pair_list(slimgb_alg* c);
302static BOOLEAN state_is(calc_state state, const int & i, const int & j, slimgb_alg* c);
303static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, slimgb_alg* c=NULL);
304static int tgb_pair_better_gen(const void* ap,const void* bp);
305static poly redTailShort(poly h, kStrategy strat);
306static poly gcd_of_terms(poly p, ring r);
307static BOOLEAN extended_product_criterion(poly p1, poly gcd1, poly p2, poly gcd2, slimgb_alg* c);
308static poly kBucketGcd(kBucket* b, ring r);
309static void multi_reduction(red_object* los, int & losl, slimgb_alg* c);
310int slim_nsize(number n, ring r);
311sorted_pair_node* quick_pop_pair(slimgb_alg* c);
312sorted_pair_node* top_pair(slimgb_alg* c);
313sorted_pair_node** add_to_basis_ideal_quotient(poly h, slimgb_alg* c, int* ip);//, BOOLEAN new_pairs=TRUE);
314sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,slimgb_alg* c);
315int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj);
316int tgb_pair_better_gen2(const void* ap,const void* bp);
317int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev);
318//static int quality(poly p, int len, slimgb_alg* c);
319/**
320   makes on each red_object in a region a single_step
321 **/
322class reduction_step{
323 public:
324  /// we assume hat all occuring red_objects have same lm, and all
325  /// occ. lm's in r[l...u] are the same, only reductor does not occur
326  virtual void reduce(red_object* r, int l, int u);
327  //int reduction_id;
328  virtual ~reduction_step();
329  slimgb_alg* c;
330  int reduction_id;
331};
332class simple_reducer:public reduction_step{
333 public:
334  poly p;
335  kBucket_pt fill_back;
336  int p_len;
337  int reducer_deg;
338  simple_reducer(poly p, int p_len,int reducer_deg, slimgb_alg* c =NULL){
339    this->p=p;
340    this->reducer_deg=reducer_deg;
341    assume(p_len==pLength(p));
342    this->p_len=p_len;
343    this->c=c;
344  }
345  virtual void pre_reduce(red_object* r, int l, int u);
346  virtual void reduce(red_object* r, int l, int u);
347  ~simple_reducer();
348
349
350  virtual void do_reduce(red_object & ro);
351};
352
353//class sum_canceling_reducer:public reduction_step {
354//  void reduce(red_object* r, int l, int u);
355//};
356struct find_erg{
357  poly expand;
358  int expand_length;
359  int to_reduce_u;
360  int to_reduce_l;
361  int reduce_by;//index of reductor
362  BOOLEAN fromS;//else from los
363
364};
365
366static void multi_reduce_step(find_erg & erg, red_object* r, slimgb_alg* c);
367static void finalize_reduction_step(reduction_step* r);
368
369template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set){
370  //Print("POSHELER:%d",sizeof(wlen_type));
371  int length=strat->sl;
372  int i;
373  int an = 0;
374  int en= length;
375
376  if ((len>setL[length])
377      || ((len==setL[length]) && (pLmCmp(set[length],p)== -1)))
378    return length+1;
379
380  loop
381  {
382    if (an >= en-1)
383    {
384      if ((len<setL[an])
385          || ((len==setL[an]) && (pLmCmp(set[an],p) == 1))) return an;
386      return en;
387    }
388    i=(an+en) / 2;
389    if ((len<setL[i])
390        || ((len==setL[i]) && (pLmCmp(set[i],p) == 1))) en=i;
391    //else if ((len>setL[i])
392    //|| ((len==setL[i]) && (pLmCmp(set[i],p) == -1))) an=i;
393    else an=i;
394  }
395
396}
397#ifdef NORO_CACHE
398#define slim_prec_cast(a) (unsigned int) (unsigned long) (a)
399#define F4mat_to_number_type(a) (number_type) slim_prec_cast(a)
400typedef unsigned short tgb_uint16;
401typedef unsigned char tgb_uint8;
402typedef unsigned int tgb_uint32;
403class NoroCacheNode{
404public:
405  NoroCacheNode** branches;
406  int branches_len;
407
408 
409  NoroCacheNode(){
410    branches=NULL;
411    branches_len=0;
412   
413  }
414  NoroCacheNode* setNode(int branch, NoroCacheNode* node){
415    if (branch>=branches_len){
416      if (branches==NULL){
417        branches_len=branch+1;
418        branches_len=si_max(branches_len,3);
419        branches=(NoroCacheNode**) omalloc(branches_len*sizeof(NoroCacheNode*));
420        int i;
421        for(i=0;i<branches_len;i++){
422          branches[i]=NULL;
423        }
424      }else{
425        int branches_len_old=branches_len;
426        branches_len=branch+1;
427        branches=(NoroCacheNode**) omrealloc(branches,branches_len*sizeof(NoroCacheNode*));
428        int i;
429        for(i=branches_len_old;i<branches_len;i++){
430          branches[i]=NULL;
431        }
432      }
433    }
434    assume(branches[branch]==NULL);
435    branches[branch]=node;
436    return node;
437  }
438  NoroCacheNode* getBranch(int branch){
439    if (branch<branches_len) return branches[branch];
440    return NULL;
441  }
442  virtual ~NoroCacheNode(){
443    int i;
444    for(i=0;i<branches_len;i++){
445      delete branches[i];
446    }
447    omfree(branches);
448  }
449  NoroCacheNode* getOrInsertBranch(int branch){
450    if ((branch<branches_len)&&(branches[branch]))
451      return branches[branch];
452    else{
453      return setNode(branch,new NoroCacheNode());
454    }
455  }
456};
457class DenseRow{
458public:
459  number* array;
460  int begin;
461  int end;
462  DenseRow(){
463    array=NULL;
464  }
465  ~DenseRow(){
466    omfree(array);
467  }
468};
469template <class number_type> class SparseRow{
470public:
471  int* idx_array;
472  number_type* coef_array;
473  int len;
474  SparseRow(){
475    len=0;
476    idx_array=NULL;
477    coef_array=NULL;
478  }
479  SparseRow<number_type>(int n){
480    len=n;
481    idx_array=(int*) omalloc(n*sizeof(int));
482    coef_array=(number_type*) omalloc(n*sizeof(number_type));
483  }
484  ~SparseRow<number_type>(){
485    omfree(idx_array);
486    omfree(coef_array);
487  }
488};
489
490template <class number_type> class DataNoroCacheNode:public NoroCacheNode{
491public:
492 
493  int value_len;
494  poly value_poly;
495  #ifdef NORO_SPARSE_ROWS_PRE
496  SparseRow<number_type>* row;
497  #else
498  DenseRow* row;
499  #endif
500  int term_index;
501  DataNoroCacheNode(poly p, int len){
502    value_len=len;
503    value_poly=p;
504    row=NULL;
505    term_index=-1;
506  }
507  #ifdef NORO_SPARSE_ROWS_PRE
508  DataNoroCacheNode(SparseRow<number_type>* row){
509    if (row!=NULL)
510      value_len=row->len;
511    else
512      value_len=0;
513    value_poly=NULL;
514    this->row=row;
515    term_index=-1;
516  }
517  #endif
518  ~DataNoroCacheNode(
519  ){
520    //p_Delete(&value_poly,currRing);
521    if (row) delete row;
522  }
523};
524template <class number_type> class TermNoroDataNode{
525public:
526  DataNoroCacheNode<number_type>* node;
527  poly t;
528};
529
530template <class number_type> class NoroCache{
531public:
532  poly temp_term;
533#ifndef NORO_NON_POLY
534  void evaluatePlaceHolder(number* row,std::vector<NoroPlaceHolder>& place_holders);
535  void evaluateRows();
536  void evaluateRows(int level, NoroCacheNode* node);
537#endif
538  void collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type>* >& res);
539  void collectIrreducibleMonomials(int level,  NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>* >& res);
540
541#ifdef NORO_RED_ARRAY_RESERVER
542  int reserved;
543  poly* recursionPolyBuffer;
544#endif
545  static const int backLinkCode=-222;
546  DataNoroCacheNode<number_type>* insert(poly term, poly nf, int len){
547    //assume(impl.find(p_Copy(term,currRing))==impl.end());
548    //assume(len==pLength(nf));
549    assume(npIsOne(p_GetCoeff(term,currRing)));
550    if (term==nf){
551      term=p_Copy(term,currRing);
552
553      ressources.push_back(term);
554      nIrreducibleMonomials++;
555      return treeInsertBackLink(term);
556     
557    } else{
558     
559      if (nf){
560        //nf=p_Copy(nf,currRing);
561        assume(p_LmCmp(nf,term,currRing)==-1);
562        ressources.push_back(nf);
563      }
564      return treeInsert(term,nf,len);
565     
566    }
567   
568    //impl[term]=std::pair<PolySimple,int> (nf,len);
569  }
570  #ifdef NORO_SPARSE_ROWS_PRE
571  DataNoroCacheNode<number_type>* insert(poly term, SparseRow<number_type>* srow){
572    //assume(impl.find(p_Copy(term,currRing))==impl.end());
573    //assume(len==pLength(nf));
574
575      return treeInsert(term,srow);
576     
577 
578    //impl[term]=std::pair<PolySimple,int> (nf,len);
579  }
580  #endif
581  DataNoroCacheNode<number_type>* insertAndTransferOwnerShip(poly t, ring r){
582   
583    ressources.push_back(t);
584    DataNoroCacheNode<number_type>* res=treeInsertBackLink(t);
585    res->term_index=nIrreducibleMonomials;
586    nIrreducibleMonomials++;
587    return res;
588  }
589  poly lookup(poly term, BOOLEAN& succ, int & len);
590  DataNoroCacheNode<number_type>* getCacheReference(poly term);
591  NoroCache(){
592    buffer=NULL;
593#ifdef NORO_RED_ARRAY_RESERVER
594    reserved=0;
595    recursionPolyBuffer=(poly*)omalloc(1000000*sizeof(poly));
596#endif
597    nIrreducibleMonomials=0;
598    nReducibleMonomials=0;
599    temp_term=pOne();
600    tempBufferSize=3000;
601    tempBuffer=omalloc(tempBufferSize);
602  }
603  void ensureTempBufferSize(size_t size){
604    if (tempBufferSize<size){
605      tempBufferSize=2*size;
606      omfree(tempBuffer);
607      tempBuffer=omalloc(tempBufferSize);
608    }
609  }
610#ifdef NORO_RED_ARRAY_RESERVER
611  poly* reserve(int n){
612    poly* res=recursionPolyBuffer+reserved;
613    reserved+=n;
614    return res;
615  }
616  void free(int n){
617    reserved-=n;
618  }
619#endif
620  ~NoroCache(){
621    int s=ressources.size();
622    int i;
623    for(i=0;i<s;i++){
624      p_Delete(&ressources[i].impl,currRing);
625    }
626    p_Delete(&temp_term,currRing);
627#ifdef NORO_RED_ARRAY_RESERVER
628    omfree(recursionPolyBuffer);
629#endif
630   omfree(tempBuffer);
631  }
632 
633  int nIrreducibleMonomials;
634  int nReducibleMonomials;
635  void* tempBuffer;
636  size_t tempBufferSize;
637protected:
638  DataNoroCacheNode<number_type>* treeInsert(poly term,poly nf,int len){
639    int i;
640    nReducibleMonomials++;
641    int nvars=pVariables;
642    NoroCacheNode* parent=&root;
643    for(i=1;i<nvars;i++){
644      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
645    }
646    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(nf,len));
647  }
648  #ifdef NORO_SPARSE_ROWS_PRE
649  DataNoroCacheNode<number_type>* treeInsert(poly term,SparseRow<number_type>* srow){
650    int i;
651    nReducibleMonomials++;
652    int nvars=pVariables;
653    NoroCacheNode* parent=&root;
654    for(i=1;i<nvars;i++){
655      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
656    }
657    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(srow));
658  }
659  #endif
660  DataNoroCacheNode<number_type>* treeInsertBackLink(poly term){
661    int i;
662    int nvars=pVariables;
663    NoroCacheNode* parent=&root;
664    for(i=1;i<nvars;i++){
665      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
666    }
667    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(term,backLinkCode));
668  }
669 
670  //@TODO descruct nodes;
671  typedef std::vector<PolySimple> poly_vec;
672  poly_vec ressources;
673  //typedef std::map<PolySimple,std::pair<PolySimple,int> > cache_map;
674  //cache_map impl;
675  NoroCacheNode root;
676  number* buffer;
677};
678template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c);
679template<class number_type> MonRedResNP<number_type> noro_red_mon_to_non_poly(poly t,  NoroCache<number_type> * cache,slimgb_alg* c)
680{
681  MonRedResNP<number_type> res_holder;
682
683
684    DataNoroCacheNode<number_type>* ref=cache->getCacheReference(t);
685    if (ref!=NULL){
686
687
688      res_holder.coef=p_GetCoeff(t,c->r);
689     
690      res_holder.ref=ref;
691      p_Delete(&t,c->r);
692      return res_holder;
693    }
694 
695  unsigned long sev=p_GetShortExpVector(t,currRing);
696  int i=kFindDivisibleByInS_easy(c->strat,t,sev);
697  if (i>=0){
698    number coef_bak=p_GetCoeff(t,c->r);
699
700    p_SetCoeff(t,npInit(1),c->r);
701    assume(npIsOne(p_GetCoeff(c->strat->S[i],c->r)));
702    number coefstrat=p_GetCoeff(c->strat->S[i],c->r);
703
704
705    poly exp_diff=cache->temp_term;
706    p_ExpVectorDiff(exp_diff,t,c->strat->S[i],c->r);
707    p_SetCoeff(exp_diff,npNeg(npInvers(coefstrat)),c->r);
708    p_Setm(exp_diff,c->r);
709    assume(c->strat->S[i]!=NULL);
710
711    poly res;
712    res=pp_Mult_mm(pNext(c->strat->S[i]),exp_diff,c->r);
713
714    int len=c->strat->lenS[i]-1;
715    SparseRow<number_type>* srow;
716    srow=noro_red_to_non_poly_t<number_type>(res,len,cache,c);
717    ref=cache->insert(t,srow);
718    p_Delete(&t,c->r);
719
720
721    res_holder.coef=coef_bak;
722    res_holder.ref=ref;
723    return res_holder;
724
725  } else {
726    number coef_bak=p_GetCoeff(t,c->r);
727    number one=npInit(1);
728    p_SetCoeff(t,one,c->r);
729 
730    res_holder.ref=cache->insertAndTransferOwnerShip(t,c->r);
731    assume(res_holder.ref!=NULL);
732    res_holder.coef=coef_bak;
733   
734    return res_holder;
735   
736  }
737
738}
739/*
740poly tree_add(poly* a,int begin, int end,ring r){
741  int d=end-begin;
742  switch(d){
743    case 0:
744      return NULL;
745    case 1:
746      return a[begin];
747    case 2:
748      return p_Add_q(a[begin],a[begin+1],r);
749    default:
750      int s=d/2;
751      return p_Add_q(tree_add(a,begin,begin+s,r),tree_add(a,begin+s,end,r),r);
752  }
753}
754*/
755#ifdef __GNUC__
756#define LIKELY(expression) (__builtin_expect(!!(expression), 1))
757#define UNLIKELY(expression) (__builtin_expect(!!(expression), 0))
758#else
759#define LIKELY(expression) (expression)
760#define UNLIKELY(expression) (expression)
761#endif
762
763template<class number_type> SparseRow<number_type>* convert_to_sparse_row(number_type* temp_array,int temp_size,int non_zeros){
764SparseRow<number_type>* res=new SparseRow<number_type>(non_zeros);
765//int pos=0;
766//Print("denseness:%f\n",((double) non_zeros/(double) temp_size));
767number_type* it_coef=res->coef_array;
768int* it_idx=res->idx_array;
769#if 0
770for(i=0;i<cache->nIrreducibleMonomials;i++){
771  if (!(0==temp_array[i])){
772 
773    res->idx_array[pos]=i;
774    res->coef_array[pos]=temp_array[i];
775
776    pos++;
777    non_zeros--;
778    if (non_zeros==0) break;
779  }
780 
781}
782#else
783int64* start=(int64*) ((void*)temp_array);
784int64* end;
785const int multiple=sizeof(int64)/sizeof(number_type);
786if (temp_size==0) end=start;
787
788else
789{ 
790  int temp_size_rounded=temp_size+(multiple-(temp_size%multiple));
791  assume(temp_size_rounded>=temp_size);
792  assume(temp_size_rounded%multiple==0);
793  assume(temp_size_rounded<temp_size+multiple);
794  number_type* nt_end=temp_array+temp_size_rounded;
795  end=(int64*)((void*)nt_end);
796}
797int64* it=start;
798while(it!=end){
799  if UNLIKELY((*it)!=0){
800    int small_i;
801    const int temp_index=((number_type*)((void*) it))-temp_array;
802    const int bound=temp_index+multiple;
803    number_type c;
804    for(small_i=temp_index;small_i<bound;small_i++){
805      if((c=temp_array[small_i])!=0){
806        //res->idx_array[pos]=small_i;
807        //res->coef_array[pos]=temp_array[small_i];
808        (*(it_idx++))=small_i;
809        (*(it_coef++))=c;
810        //pos++;
811        non_zeros--;
812       
813      }
814      if UNLIKELY(non_zeros==0) break;
815    }
816   
817  }
818  ++it;
819}
820#endif
821return res;
822}
823template <class number_type> void add_coef_times_sparse(number_type* const temp_array,int temp_size,SparseRow<number_type>* row, number coef){
824  int j;
825  number_type* const coef_array=row->coef_array;
826  int* const idx_array=row->idx_array;
827  const int len=row->len;
828  for(j=0;j<len;j++){
829    int idx=idx_array[j];
830    assume(!(npIsZero(coef)));
831    assume(!(npIsZero((number) coef_array[j])));
832    temp_array[idx]=F4mat_to_number_type(npAddM((number) temp_array[idx],npMultM((number) coef_array[j],coef)));
833    assume(idx<temp_size);
834  }
835}
836template <class number_type> void add_sparse(number_type* const temp_array,int temp_size,SparseRow<number_type>* row){
837  int j;
838
839          number_type* const coef_array=row->coef_array;
840          int* const idx_array=row->idx_array;
841          const int len=row->len;
842        for(j=0;j<len;j++){
843          int idx=idx_array[j];
844          temp_array[idx]=F4mat_to_number_type(   npAddM((number) temp_array[idx],(number) coef_array[j]));
845          assume(idx<temp_size);
846        }
847}
848template <class number_type> void sub_sparse(number_type* const temp_array,int temp_size,SparseRow<number_type>* row){
849  int j;
850
851          number_type* const coef_array=row->coef_array;
852          int* const idx_array=row->idx_array;
853          const int len=row->len;
854        for(j=0;j<len;j++){
855          int idx=idx_array[j];
856          temp_array[idx]=F4mat_to_number_type(   npSubM((number) temp_array[idx],(number) coef_array[j]));
857          assume(idx<temp_size);
858        }
859}
860template <class number_type> SparseRow<number_type>* noro_red_to_non_poly_dense(MonRedResNP<number_type>* mon, int len,NoroCache<number_type>* cache){
861  size_t temp_size_bytes=cache->nIrreducibleMonomials*sizeof(number_type)+8;//use 8bit int for testing
862   assume(sizeof(int64)==8);
863   cache->ensureTempBufferSize(temp_size_bytes);
864   number_type* temp_array=(number_type*) cache->tempBuffer;//omalloc(cache->nIrreducibleMonomials*sizeof(number_type));
865   int temp_size=cache->nIrreducibleMonomials;
866   memset(temp_array,0,temp_size_bytes);
867   number minus_one=npInit(-1);
868   int i;
869   for(i=0;i<len;i++){
870     MonRedResNP<number_type> red=mon[i];
871     if ((red.ref)){
872       if (red.ref->row){
873         SparseRow<number_type>* row=red.ref->row;
874         number coef=red.coef;
875         int j;
876         if (!((coef==(number) 1)||(coef==minus_one))){
877           add_coef_times_sparse(temp_array,temp_size,row,coef);
878
879
880
881         }else{
882           if (coef==(number) 1){
883              add_sparse(temp_array,temp_size,row);
884           } else {
885
886             sub_sparse(temp_array,temp_size,row);
887           }
888         }
889       }
890       else{
891         if (red.ref->value_len==NoroCache<number_type>::backLinkCode){
892           temp_array[red.ref->term_index]=F4mat_to_number_type( npAddM((number) temp_array[red.ref->term_index],red.coef));
893         } else {
894           //PrintS("third case\n");
895         }
896       }
897     }
898   }
899   int non_zeros=0;
900   for(i=0;i<cache->nIrreducibleMonomials;i++){
901     //if (!(temp_array[i]==0)){
902     //  non_zeros++;
903     //}
904     assume(((temp_array[i]!=0)==0)|| (((temp_array[i]!=0)==1)));
905     non_zeros+=(temp_array[i]!=0);
906   }
907
908   if (non_zeros==0){
909     //omfree(mon);
910     return NULL;
911   }
912   SparseRow<number_type>* res=convert_to_sparse_row(temp_array,temp_size, non_zeros);
913
914   //omfree(temp_array);
915
916   
917   return res;
918}
919template<class number_type> class CoefIdx{
920public:
921  number_type coef;
922  int idx;
923  bool operator<(const CoefIdx<number_type>& other) const{
924    return (idx<other.idx);
925  }
926};
927template<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){
928  int j;
929  for(j=0;j<rlen;j++){
930    assume(coef_array[j]!=0);
931    CoefIdx<number_type> ci;
932    ci.coef=F4mat_to_number_type(npMultM((number) coef,(number) coef_array[j]));
933    ci.idx=idx_array[j];
934    pairs[pos++]=ci;
935  }
936}
937template<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){
938  int j;
939  for(j=0;j<rlen;j++){
940    assume(coef_array[j]!=0);
941    CoefIdx<number_type> ci;
942    ci.coef=coef_array[j];
943    ci.idx=idx_array[j];
944    pairs[pos++]=ci;
945  }
946}
947
948template<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){
949  int j;
950  for(j=0;j<rlen;j++){
951    assume(coef_array[j]!=0);
952    CoefIdx<number_type> ci;
953    ci.coef=F4mat_to_number_type(npNegM(coef_array[j]));
954    ci.idx=idx_array[j];
955    pairs[pos++]=ci;
956  }
957}
958template <class number_type> SparseRow<number_type>* noro_red_to_non_poly_sparse(MonRedResNP<number_type>* mon, int len,NoroCache<number_type>* cache){
959  int i;
960  int together=0;
961  for(i=0;i<len;i++){
962    MonRedResNP<number_type> red=mon[i];
963    if ((red.ref) &&( red.ref->row)){
964      together+=red.ref->row->len;
965    } else {if ((red.ref) &&(red.ref->value_len==NoroCache<number_type>::backLinkCode))
966      together++;
967      }
968   
969  }
970  //PrintS("here\n");
971  if (together==0) return 0;
972  //PrintS("there\n");
973  cache->ensureTempBufferSize(together*sizeof(CoefIdx<number_type>));
974  CoefIdx<number_type>* pairs=(CoefIdx<number_type>*) cache->tempBuffer; //omalloc(together*sizeof(CoefIdx<number_type>));
975  int pos=0;
976  int j;
977  const number one=npInit(1);
978  const number minus_one=npInit(-1);
979  for(i=0;i<len;i++){
980    MonRedResNP<number_type> red=mon[i];
981    if ((red.ref) &&( red.ref->row)){
982      //together+=red.ref->row->len;
983      int* idx_array=red.ref->row->idx_array;
984      number_type* coef_array=red.ref->row->coef_array;
985      int rlen=red.ref->row->len;
986      number coef=red.coef;
987      if ((coef!=one)&&(coef!=minus_one)){
988      write_coef_times_xx_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen, coef);
989    } else
990    {
991      if (coef==one){
992        write_coef_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen);
993      } else {
994        assume(coef==minus_one);
995        write_minus_coef_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen);
996      }
997    }
998    } else {
999      if ((red.ref) &&(red.ref->value_len==NoroCache<number_type>::backLinkCode)){
1000        CoefIdx<number_type> ci;
1001        ci.coef=F4mat_to_number_type(red.coef);
1002        ci.idx=red.ref->term_index;
1003        pairs[pos++]=ci;
1004      }
1005    }
1006  }
1007  assume(pos==together);
1008  std::sort(pairs,pairs+together);
1009 
1010  int act=0;
1011 
1012  assume(pairs[0].coef!=0);
1013  for(i=1;i<together;i++){
1014    if (pairs[i].idx!=pairs[act].idx){
1015      if (pairs[act].coef!=0){
1016        act=act+1;
1017      }
1018      pairs[act]=pairs[i];
1019    } else{
1020      pairs[act].coef=F4mat_to_number_type(npAddM((number)pairs[act].coef,(number)pairs[i].coef));
1021    }
1022  }
1023 
1024  if (pairs[act].coef==0){
1025   
1026    act--;
1027  }
1028  int sparse_row_len=act+1;
1029  //Print("res len:%d",sparse_row_len);
1030  if (sparse_row_len==0) {return NULL;}
1031  SparseRow<number_type>* res=new SparseRow<number_type>(sparse_row_len);
1032  {
1033    number_type* coef_array=res->coef_array;
1034    int* idx_array=res->idx_array;
1035    for(i=0;i<sparse_row_len;i++){
1036      idx_array[i]=pairs[i].idx;
1037      coef_array[i]=pairs[i].coef;
1038    }
1039  }
1040  //omfree(pairs);
1041 
1042  return res;
1043}
1044template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c){
1045  assume(len==pLength(p));
1046  poly orig_p=p;
1047  if (p==NULL) {
1048    len=0;
1049    return NULL;
1050  }
1051 
1052  number zero=npInit(0);
1053  MonRedResNP<number_type>* mon=(MonRedResNP<number_type>*) omalloc(len*sizeof(MonRedResNP<number_type>));
1054  int i=0;
1055  double max_density=0.0;
1056  while(p){
1057
1058    poly t=p;
1059    pIter(p);
1060    pNext(t)=NULL;
1061   
1062#ifndef NDEBUG
1063    number coef_debug=p_GetCoeff(t,currRing);
1064#endif
1065    MonRedResNP<number_type> red=noro_red_mon_to_non_poly(t,cache,c);
1066    if ((red.ref) && (red.ref->row)){
1067      double act_density=(double) red.ref->row->len;
1068      act_density/=(double) cache->nIrreducibleMonomials;
1069      max_density=std::max(act_density,max_density);
1070    }
1071    mon[i]=red;
1072    i++;
1073  }
1074 
1075  assume(i==len);
1076  len=i;
1077  bool dense=true;
1078  if (max_density<0.05) dense=false;
1079  if (dense){
1080    SparseRow<number_type>* res=noro_red_to_non_poly_dense(mon,len,cache);
1081    omfree(mon);
1082    return res;
1083  } else   {
1084      SparseRow<number_type>* res=noro_red_to_non_poly_sparse(mon,len,cache);
1085      omfree(mon);
1086      return res;
1087    }
1088  //in the loop before nIrreducibleMonomials increases, so position here is important
1089 
1090}
1091#endif
1092static wlen_type pair_weighted_length(int i, int j, slimgb_alg* c);
1093wlen_type pELength(poly p, ring r);
1094int terms_sort_crit(const void* a, const void* b);
1095//void simplest_gauss_modp(number* a, int nrows,int ncols);
1096// a: a[0,0],a[0,1]....a[nrows-1,ncols-1]
1097// assume: field is Zp
1098#ifdef USE_NORO
1099
1100
1101template <class number_type > void write_poly_to_row(number_type* row, poly h, poly*terms, int tn, ring r){
1102  //poly* base=row;
1103  while(h!=NULL){
1104    //Print("h:%i\n",h);
1105    number coef=p_GetCoeff(h,r);
1106    poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit);
1107    assume(ptr_to_h!=NULL);
1108    int pos=ptr_to_h-terms;
1109    row[pos]=F4mat_to_number_type(coef);
1110    //number_type_array[base+pos]=coef;
1111    pIter(h);
1112  }
1113}
1114template <class number_type > poly row_to_poly(number_type* row, poly* terms, int tn, ring r){
1115  poly h=NULL;
1116  int j;
1117  number_type zero=0;//;npInit(0);
1118  for(j=tn-1;j>=0;j--){
1119    if (!(zero==(row[j]))){
1120      poly t=terms[j];
1121      t=p_LmInit(t,r);
1122      p_SetCoeff(t,(number) row[j],r);
1123      pNext(t)=h;
1124      h=t;
1125    }
1126   
1127  }
1128  return h;
1129}
1130template <class number_type > int modP_lastIndexRow(number_type* row,int ncols){
1131  int lastIndex;
1132  const number_type zero=0;//npInit(0);
1133  for(lastIndex=ncols-1;lastIndex>=0;lastIndex--){
1134    if (!(row[lastIndex]==zero)){
1135      return lastIndex;
1136    }
1137  }
1138  return -1;
1139}
1140template <class number_type> int term_nodes_sort_crit(const void* a, const void* b){
1141  return -pLmCmp(((TermNoroDataNode<number_type>*) a)->t,((TermNoroDataNode<number_type>*) b)->t);
1142}
1143
1144template <class number_type>class ModPMatrixBackSubstProxyOnArray;
1145template <class number_type > class ModPMatrixProxyOnArray{
1146public:
1147  friend class ModPMatrixBackSubstProxyOnArray<number_type>;
1148               
1149  int ncols,nrows;
1150  ModPMatrixProxyOnArray(number_type* array, int nrows, int ncols){
1151    this->ncols=ncols;
1152    this->nrows=nrows;
1153    rows=(number_type**) omalloc(nrows*sizeof(number_type*));
1154    startIndices=(int*)omalloc(nrows*sizeof(int));
1155    int i;
1156    for(i=0;i<nrows;i++){
1157      rows[i]=array+(i*ncols);
1158      updateStartIndex(i,-1);
1159    }
1160  }
1161  ~ModPMatrixProxyOnArray(){
1162    omfree(rows);
1163    omfree(startIndices);
1164  }
1165 
1166  void permRows(int i, int j){
1167    number_type* h=rows[i];
1168    rows[i]=rows[j];
1169    rows[j]=h;
1170    int hs=startIndices[i];
1171    startIndices[i]=startIndices[j];
1172    startIndices[j]=hs;
1173  }
1174  void multiplyRow(int row, number_type coef){
1175    int i;
1176    number_type* row_array=rows[row];
1177    for(i=startIndices[row];i<ncols;i++){
1178      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
1179    }
1180  }
1181  void reduceOtherRowsForward(int r){
1182
1183    //assume rows "under r" have bigger or equal start index
1184    number_type* row_array=rows[r];
1185    number_type zero=F4mat_to_number_type(npInit(0));
1186    int start=startIndices[r];
1187    number_type coef=row_array[start];
1188    assume(start<ncols);
1189    int other_row;
1190    assume(!(npIsZero((number) row_array[start])));
1191    if (!(npIsOne((number) coef)))
1192      multiplyRow(r,F4mat_to_number_type(npInvers((number) coef)));
1193    assume(npIsOne((number) row_array[start]));
1194    int lastIndex=modP_lastIndexRow(row_array, ncols);
1195    number minus_one=npInit(-1);
1196    for (other_row=r+1;other_row<nrows;other_row++){
1197      assume(startIndices[other_row]>=start);
1198      if (startIndices[other_row]==start){
1199        int i;
1200        number_type* other_row_array=rows[other_row];
1201        number coef2=npNeg((number) other_row_array[start]);
1202        if (coef2==minus_one){
1203          for(i=start;i<=lastIndex;i++){
1204            if (row_array[i]!=zero)
1205              other_row_array[i]=F4mat_to_number_type(npSubM((number) other_row_array[i], (number) row_array[i]));
1206          }
1207      }else {
1208          //assume(FALSE);
1209          for(i=start;i<=lastIndex;i++){
1210            if (row_array[i]!=zero)
1211            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef2,(number) row_array[i]),(number) other_row_array[i]));
1212          }
1213        }
1214        updateStartIndex(other_row,start);
1215        assume(npIsZero((number) other_row_array[start]));
1216      }
1217    }
1218  }
1219  void updateStartIndex(int row,int lower_bound){
1220    number_type* row_array=rows[row];
1221    assume((lower_bound<0)||(npIsZero((number) row_array[lower_bound])));
1222    int i;
1223    //number_type zero=npInit(0);
1224    for(i=lower_bound+1;i<ncols;i++){
1225      if (!(row_array[i]==0))
1226        break;
1227    }
1228    startIndices[row]=i;
1229  }
1230  int getStartIndex(int row){
1231    return startIndices[row];
1232  }
1233  BOOLEAN findPivot(int &r, int &c){
1234    //row>=r, col>=c
1235   
1236    while(c<ncols){
1237      int i;
1238      for(i=r;i<nrows;i++){
1239        assume(startIndices[i]>=c);
1240        if (startIndices[i]==c){
1241          //r=i;
1242          if (r!=i)
1243            permRows(r,i);
1244          return TRUE;
1245        }
1246      }
1247      c++;
1248    }
1249    return FALSE;
1250  }
1251protected:
1252  number_type** rows;
1253  int* startIndices;
1254};
1255template <class number_type > class ModPMatrixBackSubstProxyOnArray{
1256  int *startIndices;
1257  number_type** rows;
1258  int *lastReducibleIndices;
1259  int ncols;
1260  int nrows;
1261  int nonZeroUntil;
1262public:
1263  void multiplyRow(int row, number_type coef){
1264    int i;
1265    number_type* row_array=rows[row];
1266    for(i=startIndices[row];i<ncols;i++){
1267      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
1268    }
1269  }
1270  ModPMatrixBackSubstProxyOnArray<number_type> (ModPMatrixProxyOnArray<number_type> & p){
1271//  (number_type* array, int nrows, int ncols, int* startIndices, number_type** rows){
1272    //we borrow some parameters ;-)
1273    //we assume, that nobody changes the order of the rows
1274    this->startIndices=p.startIndices;
1275    this->rows=p.rows;
1276    this->ncols=p.ncols;
1277    this->nrows=p.nrows;
1278    lastReducibleIndices=(int*) omalloc(nrows*sizeof(int));
1279    nonZeroUntil=0;
1280    while(nonZeroUntil<nrows){
1281      if (startIndices[nonZeroUntil]<ncols){
1282       
1283        nonZeroUntil++;
1284      } else break;
1285     
1286    }
1287    if (TEST_OPT_PROT)
1288      Print("rank:%i\n",nonZeroUntil);
1289    nonZeroUntil--;
1290    int i;
1291    for(i=0;i<=nonZeroUntil;i++){
1292      assume(startIndices[i]<ncols);
1293      assume(!(npIsZero((number) rows[i][startIndices[i]])));
1294      assume(startIndices[i]>=i);
1295      updateLastReducibleIndex(i,nonZeroUntil+1);
1296    }
1297  }
1298  void updateLastReducibleIndex(int r, int upper_bound){
1299    number_type* row_array=rows[r];
1300    if (upper_bound>nonZeroUntil) upper_bound=nonZeroUntil+1;
1301    int i;
1302    const number_type zero=0;//npInit(0);
1303    for(i=upper_bound-1;i>r;i--){
1304      int start=startIndices[i];
1305      assume(start<ncols);
1306      if (!(row_array[start]==zero)){
1307        lastReducibleIndices[r]=start;
1308        return;
1309      }
1310    }
1311    lastReducibleIndices[r]=-1;
1312  }
1313  void backwardSubstitute(int r){
1314    int start=startIndices[r];
1315    assume(start<ncols);
1316    number_type zero=0;//npInit(0);
1317    number_type* row_array=rows[r];
1318    assume((!(npIsZero((number) row_array[start]))));
1319    assume(start<ncols);
1320    int other_row;
1321    if (!(npIsOne((number) row_array[r]))){
1322      //it should be one, but this safety is not expensive
1323      multiplyRow(r, F4mat_to_number_type(npInvers((number) row_array[start])));
1324    }
1325    int lastIndex=modP_lastIndexRow(row_array, ncols);
1326    assume(lastIndex<ncols);
1327    assume(lastIndex>=0);
1328    for(other_row=r-1;other_row>=0;other_row--){
1329      assume(lastReducibleIndices[other_row]<=start);
1330      if (lastReducibleIndices[other_row]==start){
1331        number_type* other_row_array=rows[other_row];
1332        number coef=npNeg((number) other_row_array[start]);
1333        assume(!(npIsZero(coef)));
1334        int i;
1335        assume(start>startIndices[other_row]);
1336        for(i=start;i<=lastIndex;i++){
1337          if (row_array[i]!=zero)
1338            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef,(number)row_array[i]),(number)other_row_array[i]));
1339        }
1340        updateLastReducibleIndex(other_row,r);
1341      }
1342    }
1343  }
1344  ~ModPMatrixBackSubstProxyOnArray<number_type>(){
1345    omfree(lastReducibleIndices);
1346  }
1347  void backwardSubstitute(){
1348    int i;
1349    for(i=nonZeroUntil;i>0;i--){
1350      backwardSubstitute(i);
1351    }
1352  }
1353};
1354template <class number_type > void simplest_gauss_modp(number_type* a, int nrows,int ncols){
1355  //use memmoves for changing rows
1356  if (TEST_OPT_PROT)
1357    PrintS("StartGauss\n");
1358  ModPMatrixProxyOnArray<number_type> mat(a,nrows,ncols);
1359 
1360  int c=0;
1361  int r=0;
1362  while(mat.findPivot(r,c)){
1363    //int pivot=find_pivot()
1364      mat.reduceOtherRowsForward(r);
1365    r++;
1366    c++;
1367  }
1368  ModPMatrixBackSubstProxyOnArray<number_type> backmat(mat);
1369  backmat.backwardSubstitute();
1370  //backward substitutions
1371  if (TEST_OPT_PROT)
1372    PrintS("StopGauss\n");
1373}
1374//int term_nodes_sort_crit(const void* a, const void* b);
1375template <class number_type> void noro_step(poly*p,int &pn,slimgb_alg* c){
1376  //Print("Input rows %d\n",pn);
1377  int j;
1378  if (TEST_OPT_PROT){
1379    Print("Input rows %d\n",pn);
1380  }
1381
1382  NoroCache<number_type> cache;
1383
1384  SparseRow<number_type> ** srows=(SparseRow<number_type>**) omalloc(pn*sizeof(SparseRow<number_type>*));
1385  int non_zeros=0;
1386  for(j=0;j<pn;j++){
1387   
1388    poly h=p[j];
1389    int h_len=pLength(h);
1390
1391    //number coef;
1392
1393
1394    srows[non_zeros]=noro_red_to_non_poly_t<number_type>(h,h_len,&cache,c);
1395    if (srows[non_zeros]!=NULL) non_zeros++;
1396  }
1397  std::vector<DataNoroCacheNode<number_type>*> irr_nodes;
1398  cache.collectIrreducibleMonomials(irr_nodes);
1399  //now can build up terms array
1400  //Print("historic irred Mon%d\n",cache.nIrreducibleMonomials);
1401  int n=irr_nodes.size();//cache.countIrreducibleMonomials();
1402  cache.nIrreducibleMonomials=n;
1403  if (TEST_OPT_PROT){
1404    Print("Irred Mon:%d\n",n);
1405    Print("red Mon:%d\n",cache.nReducibleMonomials);
1406  }
1407  TermNoroDataNode<number_type>* term_nodes=(TermNoroDataNode<number_type>*) omalloc(n*sizeof(TermNoroDataNode<number_type>));
1408 
1409  for(j=0;j<n;j++){
1410    assume(irr_nodes[j]!=NULL);
1411    assume(irr_nodes[j]->value_len==NoroCache<number_type>::backLinkCode);
1412    term_nodes[j].t=irr_nodes[j]->value_poly;
1413    assume(term_nodes[j].t!=NULL);
1414    term_nodes[j].node=irr_nodes[j];
1415  }
1416 
1417 
1418  qsort(term_nodes,n,sizeof(TermNoroDataNode<number_type>),term_nodes_sort_crit<number_type>);
1419  poly* terms=(poly*) omalloc(n*sizeof(poly));
1420
1421  int* old_to_new_indices=(int*) omalloc(cache.nIrreducibleMonomials*sizeof(int));
1422  for(j=0;j<n;j++){
1423    old_to_new_indices[term_nodes[j].node->term_index]=j;
1424    term_nodes[j].node->term_index=j;
1425    terms[j]=term_nodes[j].t;
1426  }
1427
1428  //if (TEST_OPT_PROT)
1429  //  Print("Evaluate Rows \n");
1430  pn=non_zeros;
1431  number_type* number_array=(number_type*) omalloc(n*pn*sizeof(number_type));
1432  memset(number_array,0,sizeof(number_type)*n*pn);
1433  number zero=npInit(0);
1434
1435  for(j=0;j<pn;j++){
1436    int i;
1437    number_type* row=number_array+n*j;
1438    /*for(i=0;i<n;i++){
1439      row[i]=zero;
1440    }*/
1441
1442    SparseRow<number_type>* srow=srows[j];
1443
1444    if (srow){
1445      int* const idx_array=srow->idx_array;
1446      number_type* const coef_array=srow->coef_array;
1447      const int len=srow->len;
1448    for(i=0;i<len;i++){
1449      int idx=old_to_new_indices[idx_array[i]];
1450      row[idx]=F4mat_to_number_type(coef_array[i]);
1451    }
1452    delete srow;
1453    }
1454  }
1455 
1456  static int export_n=0;
1457  //export_mat(number_array,pn,n,"mat%i.py",++export_n);
1458  simplest_gauss_modp(number_array,pn,n);
1459
1460  int p_pos=0;
1461  for(j=0;j<pn;j++){
1462    poly h=row_to_poly(number_array+j*n,terms,n,c->r);
1463    if(h!=NULL){
1464      p[p_pos++]=h;
1465    }
1466  }
1467  pn=p_pos;
1468  omfree(terms);
1469  omfree(term_nodes);
1470  omfree(number_array);
1471  #ifdef NORO_NON_POLY
1472  omfree(srows);
1473  omfree(old_to_new_indices);
1474  #endif
1475  //don't forget the rank
1476 
1477}
1478
1479template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type> *>& res){
1480  int i;
1481  for(i=0;i<root.branches_len;i++){
1482    collectIrreducibleMonomials(1,root.branches[i],res);
1483  }
1484}
1485template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials(int level, NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>*>& res){
1486  assume(level>=0);
1487  if (node==NULL) return;
1488  if (level<pVariables){
1489    int i,sum;
1490    for(i=0;i<node->branches_len;i++){
1491      collectIrreducibleMonomials(level+1,node->branches[i],res);
1492    }
1493  } else {
1494    DataNoroCacheNode<number_type>* dn=(DataNoroCacheNode<number_type>*) node;
1495    if (dn->value_len==backLinkCode){
1496      res.push_back(dn);
1497    } 
1498  }
1499}
1500
1501template<class number_type> DataNoroCacheNode<number_type>* NoroCache<number_type>::getCacheReference(poly term){
1502  int i;
1503  NoroCacheNode* parent=&root;
1504  for(i=1;i<pVariables;i++){
1505    parent=parent->getBranch(p_GetExp(term,i,currRing));
1506    if (!(parent)){
1507      return NULL;
1508    }
1509  }
1510  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1511  return res_holder;
1512}
1513template<class number_type> poly NoroCache<number_type>::lookup(poly term, BOOLEAN& succ, int & len){
1514  int i;
1515  NoroCacheNode* parent=&root;
1516  for(i=1;i<pVariables;i++){
1517    parent=parent->getBranch(p_GetExp(term,i,currRing));
1518    if (!(parent)){
1519      succ=FALSE;
1520      return NULL;
1521    }
1522  }
1523  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1524  if (res_holder){
1525    succ=TRUE;
1526    if ((res_holder->value_len==backLinkCode)){
1527      len=1;
1528      return term;
1529    }
1530    len=res_holder->value_len;
1531    return res_holder->value_poly;
1532  } else {
1533    succ=FALSE;
1534    return NULL;
1535  }
1536}
1537#endif
1538
1539#endif
Note: See TracBrowser for help on using the repository browser.