source: git/kernel/tgb_internal.h @ 2ade7a2

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