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

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