source: git/kernel/tgb_internal.h @ 737a68

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