source: git/kernel/tgb_internal.h @ 2320884

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