source: git/kernel/tgb_internal.h @ 9b3700

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