source: git/kernel/GBEngine/tgb_internal.h @ bd0e07

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