source: git/kernel/tgb_internal.h @ 52e2f6

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