source: git/kernel/tgb_internal.h @ 342c80

fieker-DuValspielwiese
Last change on this file since 342c80 was 86aa6a1, checked in by Michael Brickenstein <bricken@…>, 16 years ago
+ deg in exp git-svn-id: file:///usr/local/Singular/svn/trunk@10812 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 47.6 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.72 2008-06-27 12:10:58 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,int deg_pos);
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  int deg_pos;
244  BOOLEAN use_noro;
245  BOOLEAN use_noro_last_block;
246  BOOLEAN isDifficultField;
247  BOOLEAN completed;
248  BOOLEAN is_homog;
249  BOOLEAN tailReductions;
250  BOOLEAN eliminationProblem;
251  BOOLEAN F4_mode;
252  BOOLEAN nc;
253  #ifdef TGB_RESORT_PAIRS
254  BOOLEAN used_b;
255  #endif
256  unsigned long pTotaldegree(poly p){
257      pTest(p);
258      //assume(pDeg(p,r)==::pTotaldegree(p,r));
259      assume(::pTotaldegree(p,r)==p->exp[deg_pos]);
260      return p->exp[deg_pos];
261      //return ::pTotaldegree(p,this->r);
262  }
263  int pTotaldegree_full(poly p){
264    int r=0;
265    while(p){
266      int d=this->pTotaldegree(p);
267      r=si_max(r,d);
268      pIter(p);
269    }
270    return r;
271  }
272};
273class red_object{
274 public:
275  kBucket_pt bucket;
276  poly p;
277  unsigned long sev;
278  int sugar;
279  void flatten();
280  void validate();
281  wlen_type initial_quality;
282  void adjust_coefs(number c_r, number c_ac_r);
283  wlen_type guess_quality(slimgb_alg* c);
284  int clear_to_poly();
285  void canonicalize();
286};
287
288
289enum calc_state
290  {
291    UNCALCULATED,
292    HASTREP//,
293    //UNIMPORTANT,
294    //SOONTREP
295  };
296static BOOLEAN pair_cmp(sorted_pair_node* a,sorted_pair_node* b);
297template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set);
298static int add_to_reductors(slimgb_alg* c, poly h, int len, int ecart, BOOLEAN simplified=FALSE);
299static int bucket_guess(kBucket* bucket);
300static poly redNFTail (poly h,const int sl,kStrategy strat, int len);
301static poly redNF2 (poly h,slimgb_alg* c , int &len, number&  m,int n=0);
302void free_sorted_pair_node(sorted_pair_node* s, ring r);
303static void shorten_tails(slimgb_alg* c, poly monom);
304static void replace_pair(int & i, int & j, slimgb_alg* c);
305//static sorted_pair_node** add_to_basis(poly h, int i, int j,slimgb_alg* c, int* ip=NULL);
306static void do_this_spoly_stuff(int i,int j,slimgb_alg* c);
307//ideal t_rep_gb(ring r,ideal arg_I);
308ideal do_t_rep_gb(ring r,ideal arg_I, int syz_comp, BOOLEAN F4_mode,int deg_pos);
309static BOOLEAN has_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* state);
310static int* make_connections(int from, poly bound, slimgb_alg* c);
311static int* make_connections(int from, int to, poly bound, slimgb_alg* c);
312void now_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
313static void soon_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
314static int pLcmDeg(poly a, poly b);
315static int simple_posInS (kStrategy strat, poly p,int len, wlen_type wlen);
316static BOOLEAN find_next_pair(slimgb_alg* c, BOOLEAN go_higher=TRUE);
317
318static sorted_pair_node* pop_pair(slimgb_alg* c);
319static BOOLEAN no_pairs(slimgb_alg* c);
320void clean_top_of_pair_list(slimgb_alg* c);
321static void super_clean_top_of_pair_list(slimgb_alg* c);
322static BOOLEAN state_is(calc_state state, const int & i, const int & j, slimgb_alg* c);
323static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, slimgb_alg* c=NULL);
324static int tgb_pair_better_gen(const void* ap,const void* bp);
325static poly redTailShort(poly h, kStrategy strat);
326static poly gcd_of_terms(poly p, ring r);
327static BOOLEAN extended_product_criterion(poly p1, poly gcd1, poly p2, poly gcd2, slimgb_alg* c);
328static poly kBucketGcd(kBucket* b, ring r);
329static void multi_reduction(red_object* los, int & losl, slimgb_alg* c);
330int slim_nsize(number n, ring r);
331sorted_pair_node* quick_pop_pair(slimgb_alg* c);
332sorted_pair_node* top_pair(slimgb_alg* c);
333sorted_pair_node** add_to_basis_ideal_quotient(poly h, slimgb_alg* c, int* ip);//, BOOLEAN new_pairs=TRUE);
334sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,slimgb_alg* c);
335int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj);
336int tgb_pair_better_gen2(const void* ap,const void* bp);
337int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev);
338//static int quality(poly p, int len, slimgb_alg* c);
339/**
340   makes on each red_object in a region a single_step
341 **/
342class reduction_step{
343 public:
344  /// we assume hat all occuring red_objects have same lm, and all
345  /// occ. lm's in r[l...u] are the same, only reductor does not occur
346  virtual void reduce(red_object* r, int l, int u);
347  //int reduction_id;
348  virtual ~reduction_step();
349  slimgb_alg* c;
350  int reduction_id;
351};
352class simple_reducer:public reduction_step{
353 public:
354  poly p;
355  kBucket_pt fill_back;
356  int p_len;
357  int reducer_deg;
358  simple_reducer(poly p, int p_len,int reducer_deg, slimgb_alg* c =NULL){
359    this->p=p;
360    this->reducer_deg=reducer_deg;
361    assume(p_len==pLength(p));
362    this->p_len=p_len;
363    this->c=c;
364  }
365  virtual void pre_reduce(red_object* r, int l, int u);
366  virtual void reduce(red_object* r, int l, int u);
367  ~simple_reducer();
368
369
370  virtual void do_reduce(red_object & ro);
371};
372
373//class sum_canceling_reducer:public reduction_step {
374//  void reduce(red_object* r, int l, int u);
375//};
376struct find_erg{
377  poly expand;
378  int expand_length;
379  int to_reduce_u;
380  int to_reduce_l;
381  int reduce_by;//index of reductor
382  BOOLEAN fromS;//else from los
383
384};
385
386static void multi_reduce_step(find_erg & erg, red_object* r, slimgb_alg* c);
387static void finalize_reduction_step(reduction_step* r);
388
389template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set){
390  //Print("POSHELER:%d",sizeof(wlen_type));
391  int length=strat->sl;
392  int i;
393  int an = 0;
394  int en= length;
395
396  if ((len>setL[length])
397      || ((len==setL[length]) && (pLmCmp(set[length],p)== -1)))
398    return length+1;
399
400  loop
401  {
402    if (an >= en-1)
403    {
404      if ((len<setL[an])
405          || ((len==setL[an]) && (pLmCmp(set[an],p) == 1))) return an;
406      return en;
407    }
408    i=(an+en) / 2;
409    if ((len<setL[i])
410        || ((len==setL[i]) && (pLmCmp(set[i],p) == 1))) en=i;
411    //else if ((len>setL[i])
412    //|| ((len==setL[i]) && (pLmCmp(set[i],p) == -1))) an=i;
413    else an=i;
414  }
415
416}
417#ifdef NORO_CACHE
418#define slim_prec_cast(a) (unsigned int) (unsigned long) (a)
419#define F4mat_to_number_type(a) (number_type) slim_prec_cast(a)
420typedef unsigned short tgb_uint16;
421typedef unsigned char tgb_uint8;
422typedef unsigned int tgb_uint32;
423class NoroCacheNode{
424public:
425  NoroCacheNode** branches;
426  int branches_len;
427
428 
429  NoroCacheNode(){
430    branches=NULL;
431    branches_len=0;
432   
433  }
434  NoroCacheNode* setNode(int branch, NoroCacheNode* node){
435    if (branch>=branches_len){
436      if (branches==NULL){
437        branches_len=branch+1;
438        branches_len=si_max(branches_len,3);
439        branches=(NoroCacheNode**) omalloc(branches_len*sizeof(NoroCacheNode*));
440        int i;
441        for(i=0;i<branches_len;i++){
442          branches[i]=NULL;
443        }
444      }else{
445        int branches_len_old=branches_len;
446        branches_len=branch+1;
447        branches=(NoroCacheNode**) omrealloc(branches,branches_len*sizeof(NoroCacheNode*));
448        int i;
449        for(i=branches_len_old;i<branches_len;i++){
450          branches[i]=NULL;
451        }
452      }
453    }
454    assume(branches[branch]==NULL);
455    branches[branch]=node;
456    return node;
457  }
458  NoroCacheNode* getBranch(int branch){
459    if (branch<branches_len) return branches[branch];
460    return NULL;
461  }
462  virtual ~NoroCacheNode(){
463    int i;
464    for(i=0;i<branches_len;i++){
465      delete branches[i];
466    }
467    omfree(branches);
468  }
469  NoroCacheNode* getOrInsertBranch(int branch){
470    if ((branch<branches_len)&&(branches[branch]))
471      return branches[branch];
472    else{
473      return setNode(branch,new NoroCacheNode());
474    }
475  }
476};
477class DenseRow{
478public:
479  number* array;
480  int begin;
481  int end;
482  DenseRow(){
483    array=NULL;
484  }
485  ~DenseRow(){
486    omfree(array);
487  }
488};
489template <class number_type> class SparseRow{
490public:
491  int* idx_array;
492  number_type* coef_array;
493  int len;
494  SparseRow(){
495    len=0;
496    idx_array=NULL;
497    coef_array=NULL;
498  }
499  SparseRow<number_type>(int n){
500    len=n;
501    idx_array=(int*) omalloc(n*sizeof(int));
502    coef_array=(number_type*) omalloc(n*sizeof(number_type));
503  }
504  SparseRow<number_type>(int n, const number_type* source){
505    len=n;
506    idx_array=NULL;
507    coef_array=(number_type*) omalloc(n*sizeof(number_type));
508    memcpy(coef_array,source,n*sizeof(number_type));
509  }
510  ~SparseRow<number_type>(){
511    omfree(idx_array);
512    omfree(coef_array);
513  }
514};
515
516template <class number_type> class DataNoroCacheNode:public NoroCacheNode{
517public:
518 
519  int value_len;
520  poly value_poly;
521  #ifdef NORO_SPARSE_ROWS_PRE
522  SparseRow<number_type>* row;
523  #else
524  DenseRow* row;
525  #endif
526  int term_index;
527  DataNoroCacheNode(poly p, int len){
528    value_len=len;
529    value_poly=p;
530    row=NULL;
531    term_index=-1;
532  }
533  #ifdef NORO_SPARSE_ROWS_PRE
534  DataNoroCacheNode(SparseRow<number_type>* row){
535    if (row!=NULL)
536      value_len=row->len;
537    else
538      value_len=0;
539    value_poly=NULL;
540    this->row=row;
541    term_index=-1;
542  }
543  #endif
544  ~DataNoroCacheNode(
545  ){
546    //p_Delete(&value_poly,currRing);
547    if (row) delete row;
548  }
549};
550template <class number_type> class TermNoroDataNode{
551public:
552  DataNoroCacheNode<number_type>* node;
553  poly t;
554};
555
556template <class number_type> class NoroCache{
557public:
558  poly temp_term;
559#ifndef NORO_NON_POLY
560  void evaluatePlaceHolder(number* row,std::vector<NoroPlaceHolder>& place_holders);
561  void evaluateRows();
562  void evaluateRows(int level, NoroCacheNode* node);
563#endif
564  void collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type>* >& res);
565  void collectIrreducibleMonomials(int level,  NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>* >& res);
566
567#ifdef NORO_RED_ARRAY_RESERVER
568  int reserved;
569  poly* recursionPolyBuffer;
570#endif
571  static const int backLinkCode=-222;
572  DataNoroCacheNode<number_type>* insert(poly term, poly nf, int len){
573    //assume(impl.find(p_Copy(term,currRing))==impl.end());
574    //assume(len==pLength(nf));
575    assume(npIsOne(p_GetCoeff(term,currRing)));
576    if (term==nf){
577      term=p_Copy(term,currRing);
578
579      ressources.push_back(term);
580      nIrreducibleMonomials++;
581      return treeInsertBackLink(term);
582     
583    } else{
584     
585      if (nf){
586        //nf=p_Copy(nf,currRing);
587        assume(p_LmCmp(nf,term,currRing)==-1);
588        ressources.push_back(nf);
589      }
590      return treeInsert(term,nf,len);
591     
592    }
593   
594    //impl[term]=std::pair<PolySimple,int> (nf,len);
595  }
596  #ifdef NORO_SPARSE_ROWS_PRE
597  DataNoroCacheNode<number_type>* insert(poly term, SparseRow<number_type>* srow){
598    //assume(impl.find(p_Copy(term,currRing))==impl.end());
599    //assume(len==pLength(nf));
600
601      return treeInsert(term,srow);
602     
603 
604    //impl[term]=std::pair<PolySimple,int> (nf,len);
605  }
606  #endif
607  DataNoroCacheNode<number_type>* insertAndTransferOwnerShip(poly t, ring r){
608   
609    ressources.push_back(t);
610    DataNoroCacheNode<number_type>* res=treeInsertBackLink(t);
611    res->term_index=nIrreducibleMonomials;
612    nIrreducibleMonomials++;
613    return res;
614  }
615  poly lookup(poly term, BOOLEAN& succ, int & len);
616  DataNoroCacheNode<number_type>* getCacheReference(poly term);
617  NoroCache(){
618    buffer=NULL;
619#ifdef NORO_RED_ARRAY_RESERVER
620    reserved=0;
621    recursionPolyBuffer=(poly*)omalloc(1000000*sizeof(poly));
622#endif
623    nIrreducibleMonomials=0;
624    nReducibleMonomials=0;
625    temp_term=pOne();
626    tempBufferSize=3000;
627    tempBuffer=omalloc(tempBufferSize);
628  }
629  void ensureTempBufferSize(size_t size){
630    if (tempBufferSize<size){
631      tempBufferSize=2*size;
632      omfree(tempBuffer);
633      tempBuffer=omalloc(tempBufferSize);
634    }
635  }
636#ifdef NORO_RED_ARRAY_RESERVER
637  poly* reserve(int n){
638    poly* res=recursionPolyBuffer+reserved;
639    reserved+=n;
640    return res;
641  }
642  void free(int n){
643    reserved-=n;
644  }
645#endif
646  ~NoroCache(){
647    int s=ressources.size();
648    int i;
649    for(i=0;i<s;i++){
650      p_Delete(&ressources[i].impl,currRing);
651    }
652    p_Delete(&temp_term,currRing);
653#ifdef NORO_RED_ARRAY_RESERVER
654    omfree(recursionPolyBuffer);
655#endif
656   omfree(tempBuffer);
657  }
658 
659  int nIrreducibleMonomials;
660  int nReducibleMonomials;
661  void* tempBuffer;
662  size_t tempBufferSize;
663protected:
664  DataNoroCacheNode<number_type>* treeInsert(poly term,poly nf,int len){
665    int i;
666    nReducibleMonomials++;
667    int nvars=pVariables;
668    NoroCacheNode* parent=&root;
669    for(i=1;i<nvars;i++){
670      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
671    }
672    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(nf,len));
673  }
674  #ifdef NORO_SPARSE_ROWS_PRE
675  DataNoroCacheNode<number_type>* treeInsert(poly term,SparseRow<number_type>* srow){
676    int i;
677    nReducibleMonomials++;
678    int nvars=pVariables;
679    NoroCacheNode* parent=&root;
680    for(i=1;i<nvars;i++){
681      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
682    }
683    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(srow));
684  }
685  #endif
686  DataNoroCacheNode<number_type>* treeInsertBackLink(poly term){
687    int i;
688    int nvars=pVariables;
689    NoroCacheNode* parent=&root;
690    for(i=1;i<nvars;i++){
691      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
692    }
693    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(term,backLinkCode));
694  }
695 
696  //@TODO descruct nodes;
697  typedef std::vector<PolySimple> poly_vec;
698  poly_vec ressources;
699  //typedef std::map<PolySimple,std::pair<PolySimple,int> > cache_map;
700  //cache_map impl;
701  NoroCacheNode root;
702  number* buffer;
703};
704template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c);
705template<class number_type> MonRedResNP<number_type> noro_red_mon_to_non_poly(poly t,  NoroCache<number_type> * cache,slimgb_alg* c)
706{
707  MonRedResNP<number_type> res_holder;
708
709
710    DataNoroCacheNode<number_type>* ref=cache->getCacheReference(t);
711    if (ref!=NULL){
712
713
714      res_holder.coef=p_GetCoeff(t,c->r);
715     
716      res_holder.ref=ref;
717      p_Delete(&t,c->r);
718      return res_holder;
719    }
720 
721  unsigned long sev=p_GetShortExpVector(t,currRing);
722  int i=kFindDivisibleByInS_easy(c->strat,t,sev);
723  if (i>=0){
724    number coef_bak=p_GetCoeff(t,c->r);
725
726    p_SetCoeff(t,npInit(1),c->r);
727    assume(npIsOne(p_GetCoeff(c->strat->S[i],c->r)));
728    number coefstrat=p_GetCoeff(c->strat->S[i],c->r);
729
730
731    poly exp_diff=cache->temp_term;
732    p_ExpVectorDiff(exp_diff,t,c->strat->S[i],c->r);
733    p_SetCoeff(exp_diff,npNeg(npInvers(coefstrat)),c->r);
734    p_Setm(exp_diff,c->r);
735    assume(c->strat->S[i]!=NULL);
736
737    poly res;
738    res=pp_Mult_mm(pNext(c->strat->S[i]),exp_diff,c->r);
739
740    int len=c->strat->lenS[i]-1;
741    SparseRow<number_type>* srow;
742    srow=noro_red_to_non_poly_t<number_type>(res,len,cache,c);
743    ref=cache->insert(t,srow);
744    p_Delete(&t,c->r);
745
746
747    res_holder.coef=coef_bak;
748    res_holder.ref=ref;
749    return res_holder;
750
751  } else {
752    number coef_bak=p_GetCoeff(t,c->r);
753    number one=npInit(1);
754    p_SetCoeff(t,one,c->r);
755 
756    res_holder.ref=cache->insertAndTransferOwnerShip(t,c->r);
757    assume(res_holder.ref!=NULL);
758    res_holder.coef=coef_bak;
759   
760    return res_holder;
761   
762  }
763
764}
765/*
766poly tree_add(poly* a,int begin, int end,ring r){
767  int d=end-begin;
768  switch(d){
769    case 0:
770      return NULL;
771    case 1:
772      return a[begin];
773    case 2:
774      return p_Add_q(a[begin],a[begin+1],r);
775    default:
776      int s=d/2;
777      return p_Add_q(tree_add(a,begin,begin+s,r),tree_add(a,begin+s,end,r),r);
778  }
779}
780*/
781#ifdef __GNUC__
782#define LIKELY(expression) (__builtin_expect(!!(expression), 1))
783#define UNLIKELY(expression) (__builtin_expect(!!(expression), 0))
784#else
785#define LIKELY(expression) (expression)
786#define UNLIKELY(expression) (expression)
787#endif
788
789template<class number_type> SparseRow<number_type>* convert_to_sparse_row(number_type* temp_array,int temp_size,int non_zeros){
790SparseRow<number_type>* res=new SparseRow<number_type>(non_zeros);
791//int pos=0;
792//Print("denseness:%f\n",((double) non_zeros/(double) temp_size));
793number_type* it_coef=res->coef_array;
794int* it_idx=res->idx_array;
795#if 0
796for(i=0;i<cache->nIrreducibleMonomials;i++){
797  if (!(0==temp_array[i])){
798 
799    res->idx_array[pos]=i;
800    res->coef_array[pos]=temp_array[i];
801
802    pos++;
803    non_zeros--;
804    if (non_zeros==0) break;
805  }
806 
807}
808#else
809int64* start=(int64*) ((void*)temp_array);
810int64* end;
811const int multiple=sizeof(int64)/sizeof(number_type);
812if (temp_size==0) end=start;
813
814else
815{ 
816  int temp_size_rounded=temp_size+(multiple-(temp_size%multiple));
817  assume(temp_size_rounded>=temp_size);
818  assume(temp_size_rounded%multiple==0);
819  assume(temp_size_rounded<temp_size+multiple);
820  number_type* nt_end=temp_array+temp_size_rounded;
821  end=(int64*)((void*)nt_end);
822}
823int64* it=start;
824while(it!=end){
825  if UNLIKELY((*it)!=0){
826    int small_i;
827    const int temp_index=((number_type*)((void*) it))-temp_array;
828    const int bound=temp_index+multiple;
829    number_type c;
830    for(small_i=temp_index;small_i<bound;small_i++){
831      if((c=temp_array[small_i])!=0){
832        //res->idx_array[pos]=small_i;
833        //res->coef_array[pos]=temp_array[small_i];
834        (*(it_idx++))=small_i;
835        (*(it_coef++))=c;
836        //pos++;
837        non_zeros--;
838       
839      }
840      if UNLIKELY(non_zeros==0) break;
841    }
842   
843  }
844  ++it;
845}
846#endif
847return res;
848}
849template <class number_type> void add_coef_times_sparse(number_type* const temp_array,
850int temp_size,SparseRow<number_type>* row, number coef){
851  int j;
852  number_type* const coef_array=row->coef_array;
853  int* const idx_array=row->idx_array;
854  const int len=row->len;
855  tgb_uint32 buffer[256];
856  const tgb_uint32 prime=npPrimeM;
857  const tgb_uint32 c=F4mat_to_number_type(coef);
858  assume(!(npIsZero(coef)));
859  for(j=0;j<len;j=j+256){
860    const int bound=std::min(j+256,len);
861    int i;
862    int bpos=0;
863    for(i=j;i<bound;i++){
864      buffer[bpos++]=coef_array[i];
865    }
866    int bpos_bound=bound-j;
867    for(i=0;i<bpos_bound;i++){
868       buffer[i]*=c;
869     }
870    for(i=0;i<bpos_bound;i++){
871       buffer[i]=buffer[i]%prime;
872    }
873    bpos=0;
874    for(i=j;i<bound;i++){
875      int idx=idx_array[i];
876      assume(bpos<256);
877      assume(!(npIsZero((number) buffer[bpos])));
878      temp_array[idx]=F4mat_to_number_type(npAddM((number) temp_array[idx], (number) buffer[bpos++]));
879      assume(idx<temp_size);
880    }
881   
882  }
883}
884template <class number_type> void add_coef_times_dense(number_type* const temp_array,
885int temp_size,const number_type* row, int len,number coef){
886  int j;
887  const number_type* const coef_array=row;
888  //int* const idx_array=row->idx_array;
889  //const int len=temp_size;
890  tgb_uint32 buffer[256];
891  const tgb_uint32 prime=npPrimeM;
892  const tgb_uint32 c=F4mat_to_number_type(coef);
893  assume(!(npIsZero(coef)));
894  for(j=0;j<len;j=j+256){
895    const int bound=std::min(j+256,len);
896    int i;
897    int bpos=0;
898    for(i=j;i<bound;i++){
899      buffer[bpos++]=coef_array[i];
900    }
901    int bpos_bound=bound-j;
902    for(i=0;i<bpos_bound;i++){
903       buffer[i]*=c;
904     }
905    for(i=0;i<bpos_bound;i++){
906       buffer[i]=buffer[i]%prime;
907    }
908    bpos=0;
909    for(i=j;i<bound;i++){
910      //int idx=idx_array[i];
911      assume(bpos<256);
912      //assume(!(npIsZero((number) buffer[bpos])));
913      temp_array[i]=F4mat_to_number_type(npAddM((number) temp_array[i], (number) buffer[bpos++]));
914      assume(i<temp_size);
915    }
916   
917  }
918}
919template <class number_type> void add_dense(number_type* const temp_array,
920int temp_size,const number_type* row, int len){
921  //int j;
922  const number_type* const coef_array=row;
923  //int* const idx_array=row->idx_array;
924  //const int len=temp_size;
925  tgb_uint32 buffer[256];
926  const tgb_uint32 prime=npPrimeM;
927  //const tgb_uint32 c=F4mat_to_number_type(coef);
928 
929  int i;
930  for(i=0;i<len;i++){
931   
932      temp_array[i]=F4mat_to_number_type(npAddM((number) temp_array[i], (number) row[i]));
933      assume(i<temp_size);
934    }
935   
936  }
937template <class number_type> void sub_dense(number_type* const temp_array,
938int temp_size,const number_type* row, int len){
939  //int j;
940  const number_type* const coef_array=row;
941  //int* const idx_array=row->idx_array;
942  //const int len=temp_size;
943  tgb_uint32 buffer[256];
944  const tgb_uint32 prime=npPrimeM;
945  //const tgb_uint32 c=F4mat_to_number_type(coef);
946 
947  int i;
948  for(i=0;i<len;i++){
949   
950      temp_array[i]=F4mat_to_number_type(npSubM((number) temp_array[i], (number) row[i]));
951      assume(i<temp_size);
952    }
953   
954  }
955
956template <class number_type> void add_sparse(number_type* const temp_array,int temp_size,SparseRow<number_type>* row){
957  int j;
958
959          number_type* const coef_array=row->coef_array;
960          int* const idx_array=row->idx_array;
961          const int len=row->len;
962        for(j=0;j<len;j++){
963          int idx=idx_array[j];
964          temp_array[idx]=F4mat_to_number_type(   npAddM((number) temp_array[idx],(number) coef_array[j]));
965          assume(idx<temp_size);
966        }
967}
968template <class number_type> void sub_sparse(number_type* const temp_array,int temp_size,SparseRow<number_type>* row){
969  int j;
970
971          number_type* const coef_array=row->coef_array;
972          int* const idx_array=row->idx_array;
973          const int len=row->len;
974        for(j=0;j<len;j++){
975          int idx=idx_array[j];
976          temp_array[idx]=F4mat_to_number_type(   npSubM((number) temp_array[idx],(number) coef_array[j]));
977          assume(idx<temp_size);
978        }
979}
980template <class number_type> SparseRow<number_type>* noro_red_to_non_poly_dense(MonRedResNP<number_type>* mon, int len,NoroCache<number_type>* cache){
981  size_t temp_size_bytes=cache->nIrreducibleMonomials*sizeof(number_type)+8;//use 8bit int for testing
982   assume(sizeof(int64)==8);
983   cache->ensureTempBufferSize(temp_size_bytes);
984   number_type* temp_array=(number_type*) cache->tempBuffer;//omalloc(cache->nIrreducibleMonomials*sizeof(number_type));
985   int temp_size=cache->nIrreducibleMonomials;
986   memset(temp_array,0,temp_size_bytes);
987   number minus_one=npInit(-1);
988   int i;
989   for(i=0;i<len;i++){
990     MonRedResNP<number_type> red=mon[i];
991     if ((red.ref)){
992       if (red.ref->row){
993         SparseRow<number_type>* row=red.ref->row;
994         number coef=red.coef;
995         int j;
996         if (row->idx_array){
997         if (!((coef==(number) 1)||(coef==minus_one))){
998           add_coef_times_sparse(temp_array,temp_size,row,coef);
999
1000
1001
1002         }else{
1003           if (coef==(number) 1){
1004              add_sparse(temp_array,temp_size,row);
1005           } else {
1006
1007             sub_sparse(temp_array,temp_size,row);
1008           }
1009         }
1010       } else 
1011       //TODO: treat, 1,-1
1012       if (!((coef==(number) 1)||(coef==minus_one))){
1013         add_coef_times_dense(temp_array,temp_size,row->coef_array,row->len,coef);
1014       } else {
1015         if (coef==(number)1)
1016           add_dense(temp_array,temp_size,row->coef_array,row->len);
1017         else{
1018           assume(coef==minus_one);
1019           sub_dense(temp_array,temp_size,row->coef_array,row->len);
1020           //add_coef_times_dense(temp_array,temp_size,row->coef_array,row->len,coef);
1021         }
1022       }
1023       }
1024       
1025       else{
1026         if (red.ref->value_len==NoroCache<number_type>::backLinkCode){
1027           temp_array[red.ref->term_index]=F4mat_to_number_type( npAddM((number) temp_array[red.ref->term_index],red.coef));
1028         } else {
1029           //PrintS("third case\n");
1030         }
1031       }
1032     }
1033   }
1034   int non_zeros=0;
1035   for(i=0;i<cache->nIrreducibleMonomials;i++){
1036     //if (!(temp_array[i]==0)){
1037     //  non_zeros++;
1038     //}
1039     assume(((temp_array[i]!=0)==0)|| (((temp_array[i]!=0)==1)));
1040     non_zeros+=(temp_array[i]!=0);
1041   }
1042
1043   if (non_zeros==0){
1044     //omfree(mon);
1045     return NULL;
1046   }
1047   SparseRow<number_type>* res=new SparseRow<number_type>(temp_size,temp_array);//convert_to_sparse_row(temp_array,temp_size, non_zeros);
1048
1049   //omfree(temp_array);
1050
1051   
1052   return res;
1053}
1054template<class number_type> class CoefIdx{
1055public:
1056  number_type coef;
1057  int idx;
1058  bool operator<(const CoefIdx<number_type>& other) const{
1059    return (idx<other.idx);
1060  }
1061};
1062template<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){
1063  int j;
1064  for(j=0;j<rlen;j++){
1065    assume(coef_array[j]!=0);
1066    CoefIdx<number_type> ci;
1067    ci.coef=F4mat_to_number_type(npMultM((number) coef,(number) coef_array[j]));
1068    ci.idx=idx_array[j];
1069    pairs[pos++]=ci;
1070  }
1071}
1072template<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){
1073  int j;
1074 
1075  for(j=0;j<rlen;j++){
1076    if (coef_array[j]!=0){
1077    assume(coef_array[j]!=0);
1078    CoefIdx<number_type> ci;
1079    ci.coef=F4mat_to_number_type(npMultM((number) coef,(number) coef_array[j]));
1080    assume(ci.coef!=0);
1081    ci.idx=j;
1082    pairs[pos++]=ci;
1083  }
1084  }
1085}
1086template<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){
1087  int j;
1088 
1089  for(j=0;j<rlen;j++){
1090    if (coef_array[j]!=0){
1091    assume(coef_array[j]!=0);
1092    CoefIdx<number_type> ci;
1093    ci.coef=coef_array[j];
1094    assume(ci.coef!=0);
1095    ci.idx=j;
1096    pairs[pos++]=ci;
1097  }
1098  }
1099}
1100
1101template<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){
1102  int j;
1103 
1104  for(j=0;j<rlen;j++){
1105    if (coef_array[j]!=0){
1106    assume(coef_array[j]!=0);
1107    CoefIdx<number_type> ci;
1108    ci.coef=F4mat_to_number_type(npNegM((number) coef_array[j]));
1109    assume(ci.coef!=0);
1110    ci.idx=j;
1111    pairs[pos++]=ci;
1112  }
1113  }
1114}
1115template<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){
1116  int j;
1117  for(j=0;j<rlen;j++){
1118    assume(coef_array[j]!=0);
1119    CoefIdx<number_type> ci;
1120    ci.coef=coef_array[j];
1121    ci.idx=idx_array[j];
1122    pairs[pos++]=ci;
1123  }
1124}
1125
1126template<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){
1127  int j;
1128  for(j=0;j<rlen;j++){
1129    assume(coef_array[j]!=0);
1130    CoefIdx<number_type> ci;
1131    ci.coef=F4mat_to_number_type(npNegM(coef_array[j]));
1132    ci.idx=idx_array[j];
1133    pairs[pos++]=ci;
1134  }
1135}
1136template <class number_type> SparseRow<number_type>* noro_red_to_non_poly_sparse(MonRedResNP<number_type>* mon, int len,NoroCache<number_type>* cache){
1137  int i;
1138int together=0;
1139for(i=0;i<len;i++){
1140  MonRedResNP<number_type> red=mon[i];
1141  if ((red.ref) &&( red.ref->row)){
1142    together+=red.ref->row->len;
1143    } else {
1144      if ((red.ref) &&(red.ref->value_len==NoroCache<number_type>::backLinkCode))
1145      together++;
1146  }
1147
1148}
1149  //PrintS("here\n");
1150if (together==0) return 0;
1151  //PrintS("there\n");
1152cache->ensureTempBufferSize(together*sizeof(CoefIdx<number_type>));
1153CoefIdx<number_type>* pairs=(CoefIdx<number_type>*) cache->tempBuffer; //omalloc(together*sizeof(CoefIdx<number_type>));
1154int pos=0;
1155int j;
1156const number one=npInit(1);
1157const number minus_one=npInit(-1);
1158for(i=0;i<len;i++){
1159  MonRedResNP<number_type> red=mon[i];
1160  if ((red.ref) &&( red.ref->row)){
1161      //together+=red.ref->row->len;
1162    int* idx_array=red.ref->row->idx_array;
1163    number_type* coef_array=red.ref->row->coef_array;
1164    int rlen=red.ref->row->len;
1165    number coef=red.coef;
1166    if (idx_array){
1167      if ((coef!=one)&&(coef!=minus_one)){
1168        write_coef_times_xx_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen, coef);
1169        } else
1170        {
1171          if (coef==one){
1172            write_coef_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen);
1173          } else {
1174            assume(coef==minus_one);
1175            write_minus_coef_idx_to_buffer(pairs,pos,idx_array, coef_array,rlen);
1176          }
1177        }
1178      } else{
1179        if ((coef!=one)&&(coef!=minus_one)){ 
1180          write_coef_times_xx_idx_to_buffer_dense(pairs,pos,coef_array,rlen,coef);
1181        } else{
1182          if (coef==one)
1183            write_coef_idx_to_buffer_dense(pairs,pos,coef_array,rlen);
1184          else {
1185            assume(coef==minus_one);
1186            write_minus_coef_idx_to_buffer_dense(pairs,pos,coef_array,rlen);
1187          }
1188        }
1189      }
1190    }
1191    else {
1192      if ((red.ref) &&(red.ref->value_len==NoroCache<number_type>::backLinkCode)){
1193        CoefIdx<number_type> ci;
1194        ci.coef=F4mat_to_number_type(red.coef);
1195        ci.idx=red.ref->term_index;
1196        pairs[pos++]=ci;
1197      }
1198    }
1199  }
1200  assume(pos<=together);
1201  together=pos;
1202
1203  std::sort(pairs,pairs+together);
1204
1205  int act=0;
1206
1207  assume(pairs[0].coef!=0);
1208  for(i=1;i<together;i++){
1209    if (pairs[i].idx!=pairs[act].idx){
1210      if (pairs[act].coef!=0){
1211        act=act+1;
1212      }
1213      pairs[act]=pairs[i];
1214    } else{
1215      pairs[act].coef=F4mat_to_number_type(npAddM((number)pairs[act].coef,(number)pairs[i].coef));
1216    }
1217  }
1218
1219  if (pairs[act].coef==0){
1220
1221    act--;
1222  }
1223  int sparse_row_len=act+1;
1224  //Print("res len:%d",sparse_row_len);
1225  if (sparse_row_len==0) {return NULL;}
1226  SparseRow<number_type>* res=new SparseRow<number_type>(sparse_row_len);
1227  {
1228    number_type* coef_array=res->coef_array;
1229    int* idx_array=res->idx_array;
1230    for(i=0;i<sparse_row_len;i++){
1231      idx_array[i]=pairs[i].idx;
1232      coef_array[i]=pairs[i].coef;
1233    }
1234  }
1235  //omfree(pairs);
1236
1237  return res;
1238}
1239template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c){
1240  assume(len==pLength(p));
1241  poly orig_p=p;
1242  if (p==NULL) {
1243    len=0;
1244    return NULL;
1245  }
1246 
1247  number zero=npInit(0);
1248  MonRedResNP<number_type>* mon=(MonRedResNP<number_type>*) omalloc(len*sizeof(MonRedResNP<number_type>));
1249  int i=0;
1250  double max_density=0.0;
1251  while(p){
1252
1253    poly t=p;
1254    pIter(p);
1255    pNext(t)=NULL;
1256   
1257#ifndef NDEBUG
1258    number coef_debug=p_GetCoeff(t,currRing);
1259#endif
1260    MonRedResNP<number_type> red=noro_red_mon_to_non_poly(t,cache,c);
1261    if ((red.ref) && (red.ref->row)){
1262      double act_density=(double) red.ref->row->len;
1263      act_density/=(double) cache->nIrreducibleMonomials;
1264      max_density=std::max(act_density,max_density);
1265    }
1266    mon[i]=red;
1267    i++;
1268  }
1269 
1270  assume(i==len);
1271  len=i;
1272  bool dense=true;
1273  if (max_density<0.3) dense=false;
1274  if (dense){
1275    SparseRow<number_type>* res=noro_red_to_non_poly_dense(mon,len,cache);
1276    omfree(mon);
1277    return res;
1278  } else   {
1279      SparseRow<number_type>* res=noro_red_to_non_poly_sparse(mon,len,cache);
1280      omfree(mon);
1281      return res;
1282    }
1283  //in the loop before nIrreducibleMonomials increases, so position here is important
1284 
1285}
1286#endif
1287static wlen_type pair_weighted_length(int i, int j, slimgb_alg* c);
1288wlen_type pELength(poly p, ring r);
1289int terms_sort_crit(const void* a, const void* b);
1290//void simplest_gauss_modp(number* a, int nrows,int ncols);
1291// a: a[0,0],a[0,1]....a[nrows-1,ncols-1]
1292// assume: field is Zp
1293#ifdef USE_NORO
1294
1295
1296template <class number_type > void write_poly_to_row(number_type* row, poly h, poly*terms, int tn, ring r){
1297  //poly* base=row;
1298  while(h!=NULL){
1299    //Print("h:%i\n",h);
1300    number coef=p_GetCoeff(h,r);
1301    poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit);
1302    assume(ptr_to_h!=NULL);
1303    int pos=ptr_to_h-terms;
1304    row[pos]=F4mat_to_number_type(coef);
1305    //number_type_array[base+pos]=coef;
1306    pIter(h);
1307  }
1308}
1309template <class number_type > poly row_to_poly(number_type* row, poly* terms, int tn, ring r){
1310  poly h=NULL;
1311  int j;
1312  number_type zero=0;//;npInit(0);
1313  for(j=tn-1;j>=0;j--){
1314    if (!(zero==(row[j]))){
1315      poly t=terms[j];
1316      t=p_LmInit(t,r);
1317      p_SetCoeff(t,(number) row[j],r);
1318      pNext(t)=h;
1319      h=t;
1320    }
1321   
1322  }
1323  return h;
1324}
1325template <class number_type > int modP_lastIndexRow(number_type* row,int ncols){
1326  int lastIndex;
1327  const number_type zero=0;//npInit(0);
1328  for(lastIndex=ncols-1;lastIndex>=0;lastIndex--){
1329    if (!(row[lastIndex]==zero)){
1330      return lastIndex;
1331    }
1332  }
1333  return -1;
1334}
1335template <class number_type> int term_nodes_sort_crit(const void* a, const void* b){
1336  return -pLmCmp(((TermNoroDataNode<number_type>*) a)->t,((TermNoroDataNode<number_type>*) b)->t);
1337}
1338
1339template <class number_type>class ModPMatrixBackSubstProxyOnArray;
1340template <class number_type > class ModPMatrixProxyOnArray{
1341public:
1342  friend class ModPMatrixBackSubstProxyOnArray<number_type>;
1343               
1344  int ncols,nrows;
1345  ModPMatrixProxyOnArray(number_type* array, int nrows, int ncols){
1346    this->ncols=ncols;
1347    this->nrows=nrows;
1348    rows=(number_type**) omalloc(nrows*sizeof(number_type*));
1349    startIndices=(int*)omalloc(nrows*sizeof(int));
1350    int i;
1351    for(i=0;i<nrows;i++){
1352      rows[i]=array+(i*ncols);
1353      updateStartIndex(i,-1);
1354    }
1355  }
1356  ~ModPMatrixProxyOnArray(){
1357    omfree(rows);
1358    omfree(startIndices);
1359  }
1360 
1361  void permRows(int i, int j){
1362    number_type* h=rows[i];
1363    rows[i]=rows[j];
1364    rows[j]=h;
1365    int hs=startIndices[i];
1366    startIndices[i]=startIndices[j];
1367    startIndices[j]=hs;
1368  }
1369  void multiplyRow(int row, number_type coef){
1370    int i;
1371    number_type* row_array=rows[row];
1372    for(i=startIndices[row];i<ncols;i++){
1373      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
1374    }
1375  }
1376  void reduceOtherRowsForward(int r){
1377
1378    //assume rows "under r" have bigger or equal start index
1379    number_type* row_array=rows[r];
1380    number_type zero=F4mat_to_number_type(npInit(0));
1381    int start=startIndices[r];
1382    number_type coef=row_array[start];
1383    assume(start<ncols);
1384    int other_row;
1385    assume(!(npIsZero((number) row_array[start])));
1386    if (!(npIsOne((number) coef)))
1387      multiplyRow(r,F4mat_to_number_type(npInvers((number) coef)));
1388    assume(npIsOne((number) row_array[start]));
1389    int lastIndex=modP_lastIndexRow(row_array, ncols);
1390    number minus_one=npInit(-1);
1391    for (other_row=r+1;other_row<nrows;other_row++){
1392      assume(startIndices[other_row]>=start);
1393      if (startIndices[other_row]==start){
1394        int i;
1395        number_type* other_row_array=rows[other_row];
1396        number coef2=npNeg((number) other_row_array[start]);
1397        if (coef2==minus_one){
1398          for(i=start;i<=lastIndex;i++){
1399            if (row_array[i]!=zero)
1400              other_row_array[i]=F4mat_to_number_type(npSubM((number) other_row_array[i], (number) row_array[i]));
1401          }
1402      }else {
1403          //assume(FALSE);
1404          for(i=start;i<=lastIndex;i++){
1405            if (row_array[i]!=zero)
1406            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef2,(number) row_array[i]),(number) other_row_array[i]));
1407          }
1408        }
1409        updateStartIndex(other_row,start);
1410        assume(npIsZero((number) other_row_array[start]));
1411      }
1412    }
1413  }
1414  void updateStartIndex(int row,int lower_bound){
1415    number_type* row_array=rows[row];
1416    assume((lower_bound<0)||(npIsZero((number) row_array[lower_bound])));
1417    int i;
1418    //number_type zero=npInit(0);
1419    for(i=lower_bound+1;i<ncols;i++){
1420      if (!(row_array[i]==0))
1421        break;
1422    }
1423    startIndices[row]=i;
1424  }
1425  int getStartIndex(int row){
1426    return startIndices[row];
1427  }
1428  BOOLEAN findPivot(int &r, int &c){
1429    //row>=r, col>=c
1430   
1431    while(c<ncols){
1432      int i;
1433      for(i=r;i<nrows;i++){
1434        assume(startIndices[i]>=c);
1435        if (startIndices[i]==c){
1436          //r=i;
1437          if (r!=i)
1438            permRows(r,i);
1439          return TRUE;
1440        }
1441      }
1442      c++;
1443    }
1444    return FALSE;
1445  }
1446protected:
1447  number_type** rows;
1448  int* startIndices;
1449};
1450template <class number_type > class ModPMatrixBackSubstProxyOnArray{
1451  int *startIndices;
1452  number_type** rows;
1453  int *lastReducibleIndices;
1454  int ncols;
1455  int nrows;
1456  int nonZeroUntil;
1457public:
1458  void multiplyRow(int row, number_type coef){
1459    int i;
1460    number_type* row_array=rows[row];
1461    for(i=startIndices[row];i<ncols;i++){
1462      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
1463    }
1464  }
1465  ModPMatrixBackSubstProxyOnArray<number_type> (ModPMatrixProxyOnArray<number_type> & p){
1466//  (number_type* array, int nrows, int ncols, int* startIndices, number_type** rows){
1467    //we borrow some parameters ;-)
1468    //we assume, that nobody changes the order of the rows
1469    this->startIndices=p.startIndices;
1470    this->rows=p.rows;
1471    this->ncols=p.ncols;
1472    this->nrows=p.nrows;
1473    lastReducibleIndices=(int*) omalloc(nrows*sizeof(int));
1474    nonZeroUntil=0;
1475    while(nonZeroUntil<nrows){
1476      if (startIndices[nonZeroUntil]<ncols){
1477       
1478        nonZeroUntil++;
1479      } else break;
1480     
1481    }
1482    if (TEST_OPT_PROT)
1483      Print("rank:%i\n",nonZeroUntil);
1484    nonZeroUntil--;
1485    int i;
1486    for(i=0;i<=nonZeroUntil;i++){
1487      assume(startIndices[i]<ncols);
1488      assume(!(npIsZero((number) rows[i][startIndices[i]])));
1489      assume(startIndices[i]>=i);
1490      updateLastReducibleIndex(i,nonZeroUntil+1);
1491    }
1492  }
1493  void updateLastReducibleIndex(int r, int upper_bound){
1494    number_type* row_array=rows[r];
1495    if (upper_bound>nonZeroUntil) upper_bound=nonZeroUntil+1;
1496    int i;
1497    const number_type zero=0;//npInit(0);
1498    for(i=upper_bound-1;i>r;i--){
1499      int start=startIndices[i];
1500      assume(start<ncols);
1501      if (!(row_array[start]==zero)){
1502        lastReducibleIndices[r]=start;
1503        return;
1504      }
1505    }
1506    lastReducibleIndices[r]=-1;
1507  }
1508  void backwardSubstitute(int r){
1509    int start=startIndices[r];
1510    assume(start<ncols);
1511    number_type zero=0;//npInit(0);
1512    number_type* row_array=rows[r];
1513    assume((!(npIsZero((number) row_array[start]))));
1514    assume(start<ncols);
1515    int other_row;
1516    if (!(npIsOne((number) row_array[r]))){
1517      //it should be one, but this safety is not expensive
1518      multiplyRow(r, F4mat_to_number_type(npInvers((number) row_array[start])));
1519    }
1520    int lastIndex=modP_lastIndexRow(row_array, ncols);
1521    assume(lastIndex<ncols);
1522    assume(lastIndex>=0);
1523    for(other_row=r-1;other_row>=0;other_row--){
1524      assume(lastReducibleIndices[other_row]<=start);
1525      if (lastReducibleIndices[other_row]==start){
1526        number_type* other_row_array=rows[other_row];
1527        number coef=npNeg((number) other_row_array[start]);
1528        assume(!(npIsZero(coef)));
1529        int i;
1530        assume(start>startIndices[other_row]);
1531        for(i=start;i<=lastIndex;i++){
1532          if (row_array[i]!=zero)
1533            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef,(number)row_array[i]),(number)other_row_array[i]));
1534        }
1535        updateLastReducibleIndex(other_row,r);
1536      }
1537    }
1538  }
1539  ~ModPMatrixBackSubstProxyOnArray<number_type>(){
1540    omfree(lastReducibleIndices);
1541  }
1542  void backwardSubstitute(){
1543    int i;
1544    for(i=nonZeroUntil;i>0;i--){
1545      backwardSubstitute(i);
1546    }
1547  }
1548};
1549template <class number_type > void simplest_gauss_modp(number_type* a, int nrows,int ncols){
1550  //use memmoves for changing rows
1551  //if (TEST_OPT_PROT)
1552  //    PrintS("StartGauss\n");
1553  ModPMatrixProxyOnArray<number_type> mat(a,nrows,ncols);
1554 
1555  int c=0;
1556  int r=0;
1557  while(mat.findPivot(r,c)){
1558    //int pivot=find_pivot()
1559      mat.reduceOtherRowsForward(r);
1560    r++;
1561    c++;
1562  }
1563  ModPMatrixBackSubstProxyOnArray<number_type> backmat(mat);
1564  backmat.backwardSubstitute();
1565  //backward substitutions
1566  //if (TEST_OPT_PROT)
1567  //PrintS("StopGauss\n");
1568}
1569//int term_nodes_sort_crit(const void* a, const void* b);
1570template <class number_type> void noro_step(poly*p,int &pn,slimgb_alg* c){
1571  //Print("Input rows %d\n",pn);
1572  int j;
1573  if (TEST_OPT_PROT){
1574    Print("Input rows %d\n",pn);
1575  }
1576
1577  NoroCache<number_type> cache;
1578
1579  SparseRow<number_type> ** srows=(SparseRow<number_type>**) omalloc(pn*sizeof(SparseRow<number_type>*));
1580  int non_zeros=0;
1581  for(j=0;j<pn;j++){
1582   
1583    poly h=p[j];
1584    int h_len=pLength(h);
1585
1586    //number coef;
1587
1588
1589    srows[non_zeros]=noro_red_to_non_poly_t<number_type>(h,h_len,&cache,c);
1590    if (srows[non_zeros]!=NULL) non_zeros++;
1591  }
1592  std::vector<DataNoroCacheNode<number_type>*> irr_nodes;
1593  cache.collectIrreducibleMonomials(irr_nodes);
1594  //now can build up terms array
1595  //Print("historic irred Mon%d\n",cache.nIrreducibleMonomials);
1596  int n=irr_nodes.size();//cache.countIrreducibleMonomials();
1597  cache.nIrreducibleMonomials=n;
1598  if (TEST_OPT_PROT){
1599    Print("Irred Mon:%d\n",n);
1600    Print("red Mon:%d\n",cache.nReducibleMonomials);
1601  }
1602  TermNoroDataNode<number_type>* term_nodes=(TermNoroDataNode<number_type>*) omalloc(n*sizeof(TermNoroDataNode<number_type>));
1603 
1604  for(j=0;j<n;j++){
1605    assume(irr_nodes[j]!=NULL);
1606    assume(irr_nodes[j]->value_len==NoroCache<number_type>::backLinkCode);
1607    term_nodes[j].t=irr_nodes[j]->value_poly;
1608    assume(term_nodes[j].t!=NULL);
1609    term_nodes[j].node=irr_nodes[j];
1610  }
1611 
1612 
1613  qsort(term_nodes,n,sizeof(TermNoroDataNode<number_type>),term_nodes_sort_crit<number_type>);
1614  poly* terms=(poly*) omalloc(n*sizeof(poly));
1615
1616  int* old_to_new_indices=(int*) omalloc(cache.nIrreducibleMonomials*sizeof(int));
1617  for(j=0;j<n;j++){
1618    old_to_new_indices[term_nodes[j].node->term_index]=j;
1619    term_nodes[j].node->term_index=j;
1620    terms[j]=term_nodes[j].t;
1621  }
1622
1623  //if (TEST_OPT_PROT)
1624  //  Print("Evaluate Rows \n");
1625  pn=non_zeros;
1626  number_type* number_array=(number_type*) omalloc(n*pn*sizeof(number_type));
1627  memset(number_array,0,sizeof(number_type)*n*pn);
1628  number zero=npInit(0);
1629
1630  for(j=0;j<pn;j++){
1631    int i;
1632    number_type* row=number_array+n*j;
1633    /*for(i=0;i<n;i++){
1634      row[i]=zero;
1635    }*/
1636
1637    SparseRow<number_type>* srow=srows[j];
1638
1639    if (srow){
1640      int* const idx_array=srow->idx_array;
1641      number_type* const coef_array=srow->coef_array;
1642      const int len=srow->len;
1643      if (srow->idx_array){
1644        for(i=0;i<len;i++){
1645         int idx=old_to_new_indices[idx_array[i]];
1646         row[idx]=F4mat_to_number_type(coef_array[i]);
1647        } 
1648      }
1649      else {
1650        for(i=0;i<len;i++){
1651          row[old_to_new_indices[i]]=F4mat_to_number_type(coef_array[i]);
1652        }
1653      }
1654      delete srow;
1655    }
1656  }
1657 
1658  static int export_n=0;
1659  //export_mat(number_array,pn,n,"mat%i.py",++export_n);
1660  simplest_gauss_modp(number_array,pn,n);
1661
1662  int p_pos=0;
1663  for(j=0;j<pn;j++){
1664    poly h=row_to_poly(number_array+j*n,terms,n,c->r);
1665    if(h!=NULL){
1666      p[p_pos++]=h;
1667    }
1668  }
1669  pn=p_pos;
1670  omfree(terms);
1671  omfree(term_nodes);
1672  omfree(number_array);
1673  #ifdef NORO_NON_POLY
1674  omfree(srows);
1675  omfree(old_to_new_indices);
1676  #endif
1677  //don't forget the rank
1678 
1679}
1680
1681template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type> *>& res){
1682  int i;
1683  for(i=0;i<root.branches_len;i++){
1684    collectIrreducibleMonomials(1,root.branches[i],res);
1685  }
1686}
1687template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials(int level, NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>*>& res){
1688  assume(level>=0);
1689  if (node==NULL) return;
1690  if (level<pVariables){
1691    int i,sum;
1692    for(i=0;i<node->branches_len;i++){
1693      collectIrreducibleMonomials(level+1,node->branches[i],res);
1694    }
1695  } else {
1696    DataNoroCacheNode<number_type>* dn=(DataNoroCacheNode<number_type>*) node;
1697    if (dn->value_len==backLinkCode){
1698      res.push_back(dn);
1699    } 
1700  }
1701}
1702
1703template<class number_type> DataNoroCacheNode<number_type>* NoroCache<number_type>::getCacheReference(poly term){
1704  int i;
1705  NoroCacheNode* parent=&root;
1706  for(i=1;i<pVariables;i++){
1707    parent=parent->getBranch(p_GetExp(term,i,currRing));
1708    if (!(parent)){
1709      return NULL;
1710    }
1711  }
1712  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1713  return res_holder;
1714}
1715template<class number_type> poly NoroCache<number_type>::lookup(poly term, BOOLEAN& succ, int & len){
1716  int i;
1717  NoroCacheNode* parent=&root;
1718  for(i=1;i<pVariables;i++){
1719    parent=parent->getBranch(p_GetExp(term,i,currRing));
1720    if (!(parent)){
1721      succ=FALSE;
1722      return NULL;
1723    }
1724  }
1725  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1726  if (res_holder){
1727    succ=TRUE;
1728    if ((res_holder->value_len==backLinkCode)){
1729      len=1;
1730      return term;
1731    }
1732    len=res_holder->value_len;
1733    return res_holder->value_poly;
1734  } else {
1735    succ=FALSE;
1736    return NULL;
1737  }
1738}
1739#endif
1740
1741#endif
Note: See TracBrowser for help on using the repository browser.