source: git/kernel/tgb_internal.h @ 5ac8e5

spielwiese
Last change on this file since 5ac8e5 was 5ac8e5, checked in by Michael Brickenstein <bricken@…>, 17 years ago
*bricken: fully templated noro git-svn-id: file:///usr/local/Singular/svn/trunk@9889 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 35.5 KB
Line 
1#ifndef TGB_INTERNAL_H
2#define TGB_INTERNAL_H
3//!\file tgb_internal.h
4/****************************************
5*  Computer Algebra System SINGULAR     *
6****************************************/
7/* $Id: tgb_internal.h,v 1.60 2007-02-23 13:17:55 bricken Exp $ */
8/*
9 * ABSTRACT: tgb internal .h file
10*/
11#include <omalloc.h>
12#include "p_polys.h"
13
14#include "ideals.h"
15#include "ring.h"
16#include "febase.h"
17#include "structs.h"
18#include "polys.h"
19#include "stdlib.h"
20#include <modulop.h>
21//#define USE_NORO 1
22#ifdef USE_NORO
23#define NORO_CACHE 1
24#define NORO_SPARSE_ROWS_PRE 1
25#define NORO_NON_POLY 1
26#endif
27#ifdef NORO_CACHE
28//#include <map>
29#include <vector>
30#endif
31#ifdef HAVE_BOOST_DYNAMIC_BITSET_HPP
32#define  HAVE_BOOST 1
33#endif
34//#define HAVE_BOOST 1
35//#define USE_STDVECBOOL 1
36#ifdef HAVE_BOOST
37#include "boost/dynamic_bitset.hpp"
38#include <vector>
39using boost::dynamic_bitset;
40using std::vector;
41#endif
42#ifdef USE_STDVECBOOL
43#include <vector>
44using std::vector;
45#endif
46#include "kutil.h"
47#include "kInline.cc"
48#include "kstd1.h"
49#include "kbuckets.h"
50
51
52
53
54//#define TGB_DEBUG
55#define FULLREDUCTIONS
56#define HANS_IDEA
57//#define HALFREDUCTIONS
58//#define HEAD_BIN
59//#define HOMOGENEOUS_EXAMPLE
60#define REDTAIL_S
61#define PAR_N 100
62#define PAR_N_F4 5000
63#define AC_NEW_MIN 2
64#define AC_FLATTEN 1
65
66//#define FIND_DETERMINISTIC
67//#define REDTAIL_PROT
68//#define QUICK_SPOLY_TEST
69class PolySimple{
70public:
71  PolySimple(poly p){
72    impl=p;
73  }
74  PolySimple(){
75    impl=NULL;
76  }
77  PolySimple(const PolySimple& a){
78    //impl=p_Copy(a.impl,currRing);
79    impl=a.impl;
80  }
81  PolySimple& operator=(const PolySimple& p2){
82    //p_Delete(&impl,currRing);
83    //impl=p_Copy(p2.impl,currRing);
84    impl=p2.impl;
85    return *this;
86  }
87  ~PolySimple(){
88    //p_Delete(&impl,currRing);
89  }
90  bool operator< (const PolySimple& other) const{
91    return pLmCmp(impl,other.impl)<0;
92  }
93  bool operator==(const PolySimple& other){
94    return pLmEqual(impl,other.impl);
95  }
96  poly impl;
97 
98};
99template<class number_type> class DataNoroCacheNode;
100/*class MonRedRes{
101public:
102  poly p;
103  number coef;
104  BOOLEAN changed;
105  int len;
106  BOOLEAN onlyBorrowed;
107  bool operator<(const MonRedRes& other) const{
108    int cmp=p_LmCmp(p,other.p,currRing);
109    if ((cmp<0)||((cmp==0)&&((onlyBorrowed)&&(!(other.onlyBorrowed))))){
110      return true;
111    } else return false;
112  }
113  DataNoroCacheNode* ref;
114  MonRedRes(){
115    ref=NULL;
116    p=NULL;
117  }
118};*/
119template <class number_type> class MonRedResNP{
120public:
121  number coef;
122
123
124  DataNoroCacheNode<number_type>* ref;
125  MonRedResNP(){
126    ref=NULL;
127  }
128};
129struct sorted_pair_node{
130  //criterium, which is stable 0. small lcm 1. small i 2. small j
131  wlen_type expected_length;
132  poly lcm_of_lm;
133  int i;
134  int j;
135  int deg;
136 
137 
138};
139#ifdef NORO_CACHE
140#ifndef NORO_NON_POLY
141class NoroPlaceHolder{
142public:
143  DataNoroCacheNode* ref;
144  number coef;
145};
146#endif
147#endif
148//static ideal debug_Ideal;
149
150
151struct poly_list_node{
152  poly p;
153  poly_list_node* next;
154};
155
156struct int_pair_node{
157  int_pair_node* next;
158  int a;
159  int b;
160};
161struct monom_poly{
162  poly m;
163  poly f;
164};
165struct mp_array_list{
166  monom_poly* mp;
167  int size;
168  mp_array_list* next;
169};
170
171
172struct poly_array_list{
173  poly* p;
174  int size;
175  poly_array_list* next;
176};
177class slimgb_alg
178{
179  public:
180    slimgb_alg(ideal I, int syz_comp,BOOLEAN F4);
181                void introduceDelayedPairs(poly* pa,int s);
182    virtual ~slimgb_alg();
183    void cleanDegs(int lower, int upper);
184#ifndef HAVE_BOOST
185#ifdef USE_STDVECBOOL
186  vector<vector<bool> > states;
187#else
188  char** states;
189#endif
190#else
191  vector<dynamic_bitset<> > states;
192#endif
193  ideal add_later;
194  ideal S;
195  ring r;
196  int* lengths;
197  wlen_type* weighted_lengths;
198  long* short_Exps;
199  kStrategy strat;
200  int* T_deg;
201  int* T_deg_full;
202  poly tmp_lm;
203  poly* tmp_pair_lm;
204  sorted_pair_node** tmp_spn;
205  poly* expandS;
206  poly* gcd_of_terms;
207  int_pair_node* soon_free;
208  sorted_pair_node** apairs;
209  #if 0
210  BOOLEAN* modifiedS;
211  #endif
212  #ifdef TGB_RESORT_PAIRS
213  bool* replaced;
214  #endif
215  poly_list_node* to_destroy;
216  //for F4
217  mp_array_list* F;
218  poly_array_list* F_minus;
219
220  //end for F4
221#ifdef HEAD_BIN
222  struct omBin_s*   HeadBin;
223#endif
224  unsigned int reduction_steps;
225  int n;
226  //! array_lengths should be greater equal n;
227  int syz_comp;
228  int array_lengths; 
229  int normal_forms;
230  int current_degree;
231  int Rcounter;
232  int last_index;
233  int max_pairs;
234  int pair_top;
235  int easy_product_crit;
236  int extended_product_crit;
237  int average_length;
238  int lastDpBlockStart;
239  int lastCleanedDeg;
240  BOOLEAN isDifficultField;
241  BOOLEAN completed;
242  BOOLEAN is_homog;
243  BOOLEAN tailReductions;
244  BOOLEAN eliminationProblem;
245  BOOLEAN F4_mode;
246  BOOLEAN nc;
247  #ifdef TGB_RESORT_PAIRS
248  BOOLEAN used_b;
249  #endif
250};
251class red_object{
252 public:
253  kBucket_pt bucket;
254  poly p;
255  unsigned long sev;
256  int sugar;
257  void flatten();
258  void validate();
259  wlen_type initial_quality;
260  void adjust_coefs(number c_r, number c_ac_r);
261  wlen_type guess_quality(slimgb_alg* c);
262  int clear_to_poly();
263  void canonicalize();
264};
265
266
267enum calc_state
268  {
269    UNCALCULATED,
270    HASTREP//,
271    //UNIMPORTANT,
272    //SOONTREP
273  };
274static BOOLEAN pair_cmp(sorted_pair_node* a,sorted_pair_node* b);
275template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set);
276static int add_to_reductors(slimgb_alg* c, poly h, int len, int ecart, BOOLEAN simplified=FALSE);
277static int bucket_guess(kBucket* bucket);
278static poly redNFTail (poly h,const int sl,kStrategy strat, int len);
279static poly redNF2 (poly h,slimgb_alg* c , int &len, number&  m,int n=0);
280void free_sorted_pair_node(sorted_pair_node* s, ring r);
281static void shorten_tails(slimgb_alg* c, poly monom);
282static void replace_pair(int & i, int & j, slimgb_alg* c);
283//static sorted_pair_node** add_to_basis(poly h, int i, int j,slimgb_alg* c, int* ip=NULL);
284static void do_this_spoly_stuff(int i,int j,slimgb_alg* c);
285//ideal t_rep_gb(ring r,ideal arg_I);
286static BOOLEAN has_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* state);
287static int* make_connections(int from, poly bound, slimgb_alg* c);
288static int* make_connections(int from, int to, poly bound, slimgb_alg* c);
289void now_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
290static void soon_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
291static int pLcmDeg(poly a, poly b);
292static int simple_posInS (kStrategy strat, poly p,int len, wlen_type wlen);
293static BOOLEAN find_next_pair(slimgb_alg* c, BOOLEAN go_higher=TRUE);
294
295static sorted_pair_node* pop_pair(slimgb_alg* c);
296static BOOLEAN no_pairs(slimgb_alg* c);
297void clean_top_of_pair_list(slimgb_alg* c);
298static void super_clean_top_of_pair_list(slimgb_alg* c);
299static BOOLEAN state_is(calc_state state, const int & i, const int & j, slimgb_alg* c);
300static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, slimgb_alg* c=NULL);
301static int tgb_pair_better_gen(const void* ap,const void* bp);
302static poly redTailShort(poly h, kStrategy strat);
303static poly gcd_of_terms(poly p, ring r);
304static BOOLEAN extended_product_criterion(poly p1, poly gcd1, poly p2, poly gcd2, slimgb_alg* c);
305static poly kBucketGcd(kBucket* b, ring r);
306static void multi_reduction(red_object* los, int & losl, slimgb_alg* c);
307int slim_nsize(number n, ring r);
308sorted_pair_node* quick_pop_pair(slimgb_alg* c);
309sorted_pair_node* top_pair(slimgb_alg* c);
310sorted_pair_node** add_to_basis_ideal_quotient(poly h, slimgb_alg* c, int* ip);//, BOOLEAN new_pairs=TRUE);
311sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,slimgb_alg* c);
312int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj);
313int tgb_pair_better_gen2(const void* ap,const void* bp);
314int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev);
315//static int quality(poly p, int len, slimgb_alg* c);
316/**
317   makes on each red_object in a region a single_step
318 **/
319class reduction_step{
320 public:
321  /// we assume hat all occuring red_objects have same lm, and all
322  /// occ. lm's in r[l...u] are the same, only reductor does not occur
323  virtual void reduce(red_object* r, int l, int u);
324  //int reduction_id;
325  virtual ~reduction_step();
326  slimgb_alg* c;
327  int reduction_id;
328};
329class simple_reducer:public reduction_step{
330 public:
331  poly p;
332  kBucket_pt fill_back;
333  int p_len;
334  int reducer_deg;
335  simple_reducer(poly p, int p_len,int reducer_deg, slimgb_alg* c =NULL){
336    this->p=p;
337    this->reducer_deg=reducer_deg;
338    assume(p_len==pLength(p));
339    this->p_len=p_len;
340    this->c=c;
341  }
342  virtual void pre_reduce(red_object* r, int l, int u);
343  virtual void reduce(red_object* r, int l, int u);
344  ~simple_reducer();
345
346
347  virtual void do_reduce(red_object & ro);
348};
349
350//class sum_canceling_reducer:public reduction_step {
351//  void reduce(red_object* r, int l, int u);
352//};
353struct find_erg{
354  poly expand;
355  int expand_length;
356  int to_reduce_u;
357  int to_reduce_l;
358  int reduce_by;//index of reductor
359  BOOLEAN fromS;//else from los
360
361};
362
363static void multi_reduce_step(find_erg & erg, red_object* r, slimgb_alg* c);
364static void finalize_reduction_step(reduction_step* r);
365
366template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set){
367  //Print("POSHELER:%d",sizeof(wlen_type));
368  int length=strat->sl;
369  int i;
370  int an = 0;
371  int en= length;
372
373  if ((len>setL[length])
374      || ((len==setL[length]) && (pLmCmp(set[length],p)== -1)))
375    return length+1;
376
377  loop
378  {
379    if (an >= en-1)
380    {
381      if ((len<setL[an])
382          || ((len==setL[an]) && (pLmCmp(set[an],p) == 1))) return an;
383      return en;
384    }
385    i=(an+en) / 2;
386    if ((len<setL[i])
387        || ((len==setL[i]) && (pLmCmp(set[i],p) == 1))) en=i;
388    //else if ((len>setL[i])
389    //|| ((len==setL[i]) && (pLmCmp(set[i],p) == -1))) an=i;
390    else an=i;
391  }
392
393}
394#ifdef NORO_CACHE
395#define slim_prec_cast(a) (unsigned int) (a)
396#define F4mat_to_number_type(a) (number_type) slim_prec_cast(a)
397typedef unsigned short tgb_uint16;
398typedef unsigned char tgb_uint8;
399typedef unsigned int tgb_uint32;
400class NoroCacheNode{
401public:
402  NoroCacheNode** branches;
403  int branches_len;
404
405 
406  NoroCacheNode(){
407    branches=NULL;
408    branches_len=0;
409   
410  }
411  NoroCacheNode* setNode(int branch, NoroCacheNode* node){
412    if (branch>=branches_len){
413      if (branches==NULL){
414        branches_len=branch+1;
415        branches_len=si_max(branches_len,3);
416        branches=(NoroCacheNode**) omalloc(branches_len*sizeof(NoroCacheNode*));
417        int i;
418        for(i=0;i<branches_len;i++){
419          branches[i]=NULL;
420        }
421      }else{
422        int branches_len_old=branches_len;
423        branches_len=branch+1;
424        branches=(NoroCacheNode**) omrealloc(branches,branches_len*sizeof(NoroCacheNode*));
425        int i;
426        for(i=branches_len_old;i<branches_len;i++){
427          branches[i]=NULL;
428        }
429      }
430    }
431    assume(branches[branch]==NULL);
432    branches[branch]=node;
433    return node;
434  }
435  NoroCacheNode* getBranch(int branch){
436    if (branch<branches_len) return branches[branch];
437    return NULL;
438  }
439  virtual ~NoroCacheNode(){
440    int i;
441    for(i=0;i<branches_len;i++){
442      delete branches[i];
443    }
444    omfree(branches);
445  }
446  NoroCacheNode* getOrInsertBranch(int branch){
447    if ((branch<branches_len)&&(branches[branch]))
448      return branches[branch];
449    else{
450      return setNode(branch,new NoroCacheNode());
451    }
452  }
453};
454class DenseRow{
455public:
456  number* array;
457  int begin;
458  int end;
459  DenseRow(){
460    array=NULL;
461  }
462  ~DenseRow(){
463    omfree(array);
464  }
465};
466template <class number_type> class SparseRow{
467public:
468  int* idx_array;
469  number_type* coef_array;
470  int len;
471  SparseRow(){
472    len=0;
473    idx_array=NULL;
474    coef_array=NULL;
475  }
476  SparseRow<number_type>(int n){
477    len=n;
478    idx_array=(int*) omalloc(n*sizeof(int));
479    coef_array=(number_type*) omalloc(n*sizeof(number_type));
480  }
481  ~SparseRow<number_type>(){
482    omfree(idx_array);
483    omfree(coef_array);
484  }
485};
486
487template <class number_type> class DataNoroCacheNode:public NoroCacheNode{
488public:
489 
490  int value_len;
491  poly value_poly;
492  #ifdef NORO_SPARSE_ROWS_PRE
493  SparseRow<number_type>* row;
494  #else
495  DenseRow* row;
496  #endif
497  int term_index;
498  DataNoroCacheNode(poly p, int len){
499    value_len=len;
500    value_poly=p;
501    row=NULL;
502    term_index=-1;
503  }
504  #ifdef NORO_SPARSE_ROWS_PRE
505  DataNoroCacheNode(SparseRow<number_type>* row){
506    if (row!=NULL)
507      value_len=row->len;
508    else
509      value_len=0;
510    value_poly=NULL;
511    this->row=row;
512    term_index=-1;
513  }
514  #endif
515  ~DataNoroCacheNode(
516  ){
517    //p_Delete(&value_poly,currRing);
518    if (row) delete row;
519  }
520};
521template <class number_type> class TermNoroDataNode{
522public:
523  DataNoroCacheNode<number_type>* node;
524  poly t;
525};
526
527template <class number_type> class NoroCache{
528public:
529  poly temp_term;
530#ifndef NORO_NON_POLY
531  void evaluatePlaceHolder(number* row,std::vector<NoroPlaceHolder>& place_holders);
532  void evaluateRows();
533  void evaluateRows(int level, NoroCacheNode* node);
534#endif
535  void collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type>* >& res);
536  void collectIrreducibleMonomials(int level,  NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>* >& res);
537
538#ifdef NORO_RED_ARRAY_RESERVER
539  int reserved;
540  poly* recursionPolyBuffer;
541#endif
542  static const int backLinkCode=-222;
543  DataNoroCacheNode<number_type>* insert(poly term, poly nf, int len){
544    //assume(impl.find(p_Copy(term,currRing))==impl.end());
545    //assume(len==pLength(nf));
546    assume(npIsOne(p_GetCoeff(term,currRing)));
547    if (term==nf){
548      term=p_Copy(term,currRing);
549
550      ressources.push_back(term);
551      nIrreducibleMonomials++;
552      return treeInsertBackLink(term);
553     
554    } else{
555     
556      if (nf){
557        //nf=p_Copy(nf,currRing);
558        assume(p_LmCmp(nf,term,currRing)==-1);
559        ressources.push_back(nf);
560      }
561      return treeInsert(term,nf,len);
562     
563    }
564   
565    //impl[term]=std::pair<PolySimple,int> (nf,len);
566  }
567  #ifdef NORO_SPARSE_ROWS_PRE
568  DataNoroCacheNode<number_type>* insert(poly term, SparseRow<number_type>* srow){
569    //assume(impl.find(p_Copy(term,currRing))==impl.end());
570    //assume(len==pLength(nf));
571
572      return treeInsert(term,srow);
573     
574 
575    //impl[term]=std::pair<PolySimple,int> (nf,len);
576  }
577  #endif
578  DataNoroCacheNode<number_type>* insertAndTransferOwnerShip(poly t, ring r){
579   
580    ressources.push_back(t);
581    DataNoroCacheNode<number_type>* res=treeInsertBackLink(t);
582    res->term_index=nIrreducibleMonomials;
583    nIrreducibleMonomials++;
584    return res;
585  }
586  poly lookup(poly term, BOOLEAN& succ, int & len);
587  DataNoroCacheNode<number_type>* getCacheReference(poly term);
588  NoroCache(){
589    buffer=NULL;
590#ifdef NORO_RED_ARRAY_RESERVER
591    reserved=0;
592    recursionPolyBuffer=(poly*)omalloc(1000000*sizeof(poly));
593#endif
594    nIrreducibleMonomials=0;
595    nReducibleMonomials=0;
596    temp_term=pOne();
597    tempBufferSize=3000;
598    tempBuffer=omalloc(tempBufferSize);
599  }
600  void ensureTempBufferSize(size_t size){
601    if (tempBufferSize<size){
602      tempBufferSize=2*size;
603      omfree(tempBuffer);
604      tempBuffer=omalloc(tempBufferSize);
605    }
606  }
607#ifdef NORO_RED_ARRAY_RESERVER
608  poly* reserve(int n){
609    poly* res=recursionPolyBuffer+reserved;
610    reserved+=n;
611    return res;
612  }
613  void free(int n){
614    reserved-=n;
615  }
616#endif
617  ~NoroCache(){
618    int s=ressources.size();
619    int i;
620    for(i=0;i<s;i++){
621      p_Delete(&ressources[i].impl,currRing);
622    }
623    p_Delete(&temp_term,currRing);
624#ifdef NORO_RED_ARRAY_RESERVER
625    omfree(recursionPolyBuffer);
626#endif
627   omfree(tempBuffer);
628  }
629 
630  int nIrreducibleMonomials;
631  int nReducibleMonomials;
632  void* tempBuffer;
633  size_t tempBufferSize;
634protected:
635  DataNoroCacheNode<number_type>* treeInsert(poly term,poly nf,int len){
636    int i;
637    nReducibleMonomials++;
638    int nvars=pVariables;
639    NoroCacheNode* parent=&root;
640    for(i=1;i<nvars;i++){
641      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
642    }
643    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(nf,len));
644  }
645  #ifdef NORO_SPARSE_ROWS_PRE
646  DataNoroCacheNode<number_type>* treeInsert(poly term,SparseRow<number_type>* srow){
647    int i;
648    nReducibleMonomials++;
649    int nvars=pVariables;
650    NoroCacheNode* parent=&root;
651    for(i=1;i<nvars;i++){
652      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
653    }
654    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(srow));
655  }
656  #endif
657  DataNoroCacheNode<number_type>* treeInsertBackLink(poly term){
658    int i;
659    int nvars=pVariables;
660    NoroCacheNode* parent=&root;
661    for(i=1;i<nvars;i++){
662      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
663    }
664    return (DataNoroCacheNode<number_type>*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode<number_type>(term,backLinkCode));
665  }
666 
667  //@TODO descruct nodes;
668  typedef std::vector<PolySimple> poly_vec;
669  poly_vec ressources;
670  //typedef std::map<PolySimple,std::pair<PolySimple,int> > cache_map;
671  //cache_map impl;
672  NoroCacheNode root;
673  number* buffer;
674};
675template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c);
676template<class number_type> MonRedResNP<number_type> noro_red_mon_to_non_poly(poly t,  NoroCache<number_type> * cache,slimgb_alg* c)
677{
678  MonRedResNP<number_type> res_holder;
679
680
681    DataNoroCacheNode<number_type>* ref=cache->getCacheReference(t);
682    if (ref!=NULL){
683
684
685      res_holder.coef=p_GetCoeff(t,c->r);
686     
687      res_holder.ref=ref;
688      p_Delete(&t,c->r);
689      return res_holder;
690    }
691 
692  unsigned long sev=p_GetShortExpVector(t,currRing);
693  int i=kFindDivisibleByInS_easy(c->strat,t,sev);
694  if (i>=0){
695    number coef_bak=p_GetCoeff(t,c->r);
696
697    p_SetCoeff(t,npInit(1),c->r);
698    assume(npIsOne(p_GetCoeff(c->strat->S[i],c->r)));
699    number coefstrat=p_GetCoeff(c->strat->S[i],c->r);
700
701
702    poly exp_diff=cache->temp_term;
703    p_ExpVectorDiff(exp_diff,t,c->strat->S[i],c->r);
704    p_SetCoeff(exp_diff,npNeg(npInvers(coefstrat)),c->r);
705    p_Setm(exp_diff,c->r);
706    assume(c->strat->S[i]!=NULL);
707
708    poly res;
709    res=pp_Mult_mm(pNext(c->strat->S[i]),exp_diff,c->r);
710
711    int len=c->strat->lenS[i]-1;
712    SparseRow<number_type>* srow;
713    srow=noro_red_to_non_poly_t<number_type>(res,len,cache,c);
714    ref=cache->insert(t,srow);
715    p_Delete(&t,c->r);
716
717
718    res_holder.coef=coef_bak;
719    res_holder.ref=ref;
720    return res_holder;
721
722  } else {
723    number coef_bak=p_GetCoeff(t,c->r);
724    number one=npInit(1);
725    p_SetCoeff(t,one,c->r);
726 
727    res_holder.ref=cache->insertAndTransferOwnerShip(t,c->r);
728    assume(res_holder.ref!=NULL);
729    res_holder.coef=coef_bak;
730   
731    return res_holder;
732   
733  }
734
735}
736/*
737poly tree_add(poly* a,int begin, int end,ring r){
738  int d=end-begin;
739  switch(d){
740    case 0:
741      return NULL;
742    case 1:
743      return a[begin];
744    case 2:
745      return p_Add_q(a[begin],a[begin+1],r);
746    default:
747      int s=d/2;
748      return p_Add_q(tree_add(a,begin,begin+s,r),tree_add(a,begin+s,end,r),r);
749  }
750}
751*/
752template<class number_type> SparseRow<number_type> * noro_red_to_non_poly_t(poly p, int &len, NoroCache<number_type>* cache,slimgb_alg* c){
753  assume(len==pLength(p));
754  poly orig_p=p;
755  if (p==NULL) {
756    len=0;
757    return NULL;
758  }
759 
760  number zero=npInit(0);
761  MonRedResNP<number_type>* mon=(MonRedResNP<number_type>*) omalloc(len*sizeof(MonRedResNP<number_type>));
762  int i=0;
763
764  while(p){
765
766    poly t=p;
767    pIter(p);
768    pNext(t)=NULL;
769   
770#ifndef NDEBUG
771    number coef_debug=p_GetCoeff(t,currRing);
772#endif
773    MonRedResNP<number_type> red=noro_red_mon_to_non_poly(t,cache,c);
774    mon[i]=red;
775    i++;
776  }
777 
778  assume(i==len);
779  len=i;
780  //in the loop before nIrreducibleMonomials increases, so position here is important
781  size_t temp_size_bytes=cache->nIrreducibleMonomials*sizeof(number_type)+8;//use 8bit int for testing
782  assume(sizeof(int64)==8);
783  cache->ensureTempBufferSize(temp_size_bytes);
784  number_type* temp_array=(number_type*) cache->tempBuffer;//omalloc(cache->nIrreducibleMonomials*sizeof(number_type));
785  int temp_size=cache->nIrreducibleMonomials;
786  memset(temp_array,0,temp_size_bytes);
787  for(i=0;i<len;i++){
788    MonRedResNP<number_type> red=mon[i];
789    if ((red.ref)){
790      if (red.ref->row){
791        SparseRow<number_type>* row=red.ref->row;
792        number coef=red.coef;
793        int j;
794        if (!(npIsOne(coef))){
795        for(j=0;j<row->len;j++){
796          int idx=row->idx_array[j];
797          assume(!(npIsZero(coef)));
798          assume(!(npIsZero((number) row->coef_array[j])));
799          temp_array[idx]=F4mat_to_number_type(npAddM((number) temp_array[idx],npMultM((number) row->coef_array[j],coef)));
800          assume(idx<temp_size);
801        }}else{
802          for(j=0;j<row->len;j++){
803            int idx=row->idx_array[j];
804            temp_array[idx]=F4mat_to_number_type(   npAddM((number) temp_array[idx],(number) row->coef_array[j]));
805            assume(idx<temp_size);
806          }
807        }
808      }
809      else{
810        if (red.ref->value_len==NoroCache<number_type>::backLinkCode){
811          temp_array[red.ref->term_index]=(number_type) (unsigned int) npAddM((number) temp_array[red.ref->term_index],red.coef);
812        } else {
813          //PrintS("third case\n");
814        }
815      }
816    }
817  }
818  int non_zeros=0;
819  for(i=0;i<cache->nIrreducibleMonomials;i++){
820    //if (!(temp_array[i]==0)){
821    //  non_zeros++;
822    //}
823    assume(((temp_array[i]!=0)==0)|| (((temp_array[i]!=0)==1)));
824    non_zeros+=(temp_array[i]!=0);
825  }
826 
827  if (non_zeros==0){
828    omfree(mon);
829    return NULL;
830  }
831  SparseRow<number_type>* res=new SparseRow<number_type>(non_zeros);
832  int pos=0;
833  #if 0
834  for(i=0;i<cache->nIrreducibleMonomials;i++){
835    if (!(0==temp_array[i])){
836   
837      res->idx_array[pos]=i;
838      res->coef_array[pos]=temp_array[i];
839
840      pos++;
841      non_zeros--;
842      if (non_zeros==0) break;
843    }
844   
845  }
846  #else
847  int64* start=(int64*) ((void*)temp_array);
848  int64* end;
849  const int multiple=sizeof(int64)/sizeof(number_type);
850  if (temp_size==0) end=start;
851 
852  else
853  { 
854    int temp_size_rounded=temp_size+(multiple-(temp_size%multiple));
855    assume(temp_size_rounded>=temp_size);
856    assume(temp_size_rounded%multiple==0);
857    assume(temp_size_rounded<temp_size+multiple);
858    number_type* nt_end=temp_array+temp_size_rounded;
859    end=(int64*)((void*)nt_end);
860  }
861  int64* it=start;
862  while(it!=end){
863    if ((*it)!=0){
864      int small_i;
865      const int temp_index=((number_type*)((void*) it))-temp_array;
866      const int bound=temp_index+multiple;
867      for(small_i=temp_index;small_i<bound;small_i++){
868        if(temp_array[small_i]!=0){
869          res->idx_array[pos]=small_i;
870          res->coef_array[pos]=temp_array[small_i];
871
872          pos++;
873          non_zeros--;
874         
875        }
876        if (non_zeros==0) break;
877      }
878     
879    }
880    ++it;
881  }
882  #endif
883  //omfree(temp_array);
884
885  omfree(mon);
886  return res;
887}
888#endif
889static wlen_type pair_weighted_length(int i, int j, slimgb_alg* c);
890wlen_type pELength(poly p, ring r);
891int terms_sort_crit(const void* a, const void* b);
892//void simplest_gauss_modp(number* a, int nrows,int ncols);
893// a: a[0,0],a[0,1]....a[nrows-1,ncols-1]
894// assume: field is Zp
895#ifdef USE_NORO
896
897
898template <class number_type > void write_poly_to_row(number_type* row, poly h, poly*terms, int tn, ring r){
899  //poly* base=row;
900  while(h!=NULL){
901    //Print("h:%i\n",h);
902    number coef=p_GetCoeff(h,r);
903    poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit);
904    assume(ptr_to_h!=NULL);
905    int pos=ptr_to_h-terms;
906    row[pos]=F4mat_to_number_type(coef);
907    //number_type_array[base+pos]=coef;
908    pIter(h);
909  }
910}
911template <class number_type > poly row_to_poly(number_type* row, poly* terms, int tn, ring r){
912  poly h=NULL;
913  int j;
914  number_type zero=0;//;npInit(0);
915  for(j=tn-1;j>=0;j--){
916    if (!(zero==(row[j]))){
917      poly t=terms[j];
918      t=p_LmInit(t,r);
919      p_SetCoeff(t,(number) row[j],r);
920      pNext(t)=h;
921      h=t;
922    }
923   
924  }
925  return h;
926}
927template <class number_type > int modP_lastIndexRow(number_type* row,int ncols){
928  int lastIndex;
929  const number_type zero=0;//npInit(0);
930  for(lastIndex=ncols-1;lastIndex>=0;lastIndex--){
931    if (!(row[lastIndex]==zero)){
932      return lastIndex;
933    }
934  }
935  return -1;
936}
937template <class number_type> int term_nodes_sort_crit(const void* a, const void* b){
938  return -pLmCmp(((TermNoroDataNode<number_type>*) a)->t,((TermNoroDataNode<number_type>*) b)->t);
939}
940
941template <class number_type>class ModPMatrixBackSubstProxyOnArray;
942template <class number_type > class ModPMatrixProxyOnArray{
943public:
944  friend class ModPMatrixBackSubstProxyOnArray<number_type>;
945               
946  int ncols,nrows;
947  ModPMatrixProxyOnArray(number_type* array, int nrows, int ncols){
948    this->ncols=ncols;
949    this->nrows=nrows;
950    rows=(number_type**) omalloc(nrows*sizeof(number_type*));
951    startIndices=(int*)omalloc(nrows*sizeof(int));
952    int i;
953    for(i=0;i<nrows;i++){
954      rows[i]=array+(i*ncols);
955      updateStartIndex(i,-1);
956    }
957  }
958  ~ModPMatrixProxyOnArray(){
959    omfree(rows);
960    omfree(startIndices);
961  }
962 
963  void permRows(int i, int j){
964    number_type* h=rows[i];
965    rows[i]=rows[j];
966    rows[j]=h;
967    int hs=startIndices[i];
968    startIndices[i]=startIndices[j];
969    startIndices[j]=hs;
970  }
971  void multiplyRow(int row, number_type coef){
972    int i;
973    number_type* row_array=rows[row];
974    for(i=startIndices[row];i<ncols;i++){
975      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
976    }
977  }
978  void reduceOtherRowsForward(int r){
979
980    //assume rows "under r" have bigger or equal start index
981    number_type* row_array=rows[r];
982    number_type zero=F4mat_to_number_type(npInit(0));
983    int start=startIndices[r];
984    number_type coef=row_array[start];
985    assume(start<ncols);
986    int other_row;
987    assume(!(npIsZero((number) row_array[start])));
988    if (!(npIsOne((number) coef)))
989      multiplyRow(r,F4mat_to_number_type(npInvers((number) coef)));
990    assume(npIsOne((number) row_array[start]));
991    int lastIndex=modP_lastIndexRow(row_array, ncols);
992    number minus_one=npInit(-1);
993    for (other_row=r+1;other_row<nrows;other_row++){
994      assume(startIndices[other_row]>=start);
995      if (startIndices[other_row]==start){
996        int i;
997        number_type* other_row_array=rows[other_row];
998        number coef2=npNeg((number) other_row_array[start]);
999        if (coef2==minus_one){
1000          for(i=start;i<=lastIndex;i++){
1001            if (row_array[i]!=zero)
1002              other_row_array[i]=F4mat_to_number_type(npSubM((number) other_row_array[i], (number) row_array[i]));
1003          }
1004      }else {
1005          //assume(FALSE);
1006          for(i=start;i<=lastIndex;i++){
1007            if (row_array[i]!=zero)
1008            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef2,(number) row_array[i]),(number) other_row_array[i]));
1009          }
1010        }
1011        updateStartIndex(other_row,start);
1012        assume(npIsZero((number) other_row_array[start]));
1013      }
1014    }
1015  }
1016  void updateStartIndex(int row,int lower_bound){
1017    number_type* row_array=rows[row];
1018    assume((lower_bound<0)||(npIsZero((number) row_array[lower_bound])));
1019    int i;
1020    //number_type zero=npInit(0);
1021    for(i=lower_bound+1;i<ncols;i++){
1022      if (!(row_array[i]==0))
1023        break;
1024    }
1025    startIndices[row]=i;
1026  }
1027  int getStartIndex(int row){
1028    return startIndices[row];
1029  }
1030  BOOLEAN findPivot(int &r, int &c){
1031    //row>=r, col>=c
1032   
1033    while(c<ncols){
1034      int i;
1035      for(i=r;i<nrows;i++){
1036        assume(startIndices[i]>=c);
1037        if (startIndices[i]==c){
1038          //r=i;
1039          if (r!=i)
1040            permRows(r,i);
1041          return TRUE;
1042        }
1043      }
1044      c++;
1045    }
1046    return FALSE;
1047  }
1048protected:
1049  number_type** rows;
1050  int* startIndices;
1051};
1052template <class number_type > class ModPMatrixBackSubstProxyOnArray{
1053  int *startIndices;
1054  number_type** rows;
1055  int *lastReducibleIndices;
1056  int ncols;
1057  int nrows;
1058  int nonZeroUntil;
1059public:
1060  void multiplyRow(int row, number_type coef){
1061    int i;
1062    number_type* row_array=rows[row];
1063    for(i=startIndices[row];i<ncols;i++){
1064      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
1065    }
1066  }
1067  ModPMatrixBackSubstProxyOnArray<number_type> (ModPMatrixProxyOnArray<number_type> & p){
1068//  (number_type* array, int nrows, int ncols, int* startIndices, number_type** rows){
1069    //we borrow some parameters ;-)
1070    //we assume, that nobody changes the order of the rows
1071    this->startIndices=p.startIndices;
1072    this->rows=p.rows;
1073    this->ncols=p.ncols;
1074    this->nrows=p.nrows;
1075    lastReducibleIndices=(int*) omalloc(nrows*sizeof(int));
1076    nonZeroUntil=0;
1077    while(nonZeroUntil<nrows){
1078      if (startIndices[nonZeroUntil]<ncols){
1079       
1080        nonZeroUntil++;
1081      } else break;
1082     
1083    }
1084    if (TEST_OPT_PROT)
1085      Print("rank:%i\n",nonZeroUntil);
1086    nonZeroUntil--;
1087    int i;
1088    for(i=0;i<=nonZeroUntil;i++){
1089      assume(startIndices[i]<ncols);
1090      assume(!(npIsZero((number) rows[i][startIndices[i]])));
1091      assume(startIndices[i]>=i);
1092      updateLastReducibleIndex(i,nonZeroUntil+1);
1093    }
1094  }
1095  void updateLastReducibleIndex(int r, int upper_bound){
1096    number_type* row_array=rows[r];
1097    if (upper_bound>nonZeroUntil) upper_bound=nonZeroUntil+1;
1098    int i;
1099    const number_type zero=0;//npInit(0);
1100    for(i=upper_bound-1;i>r;i--){
1101      int start=startIndices[i];
1102      assume(start<ncols);
1103      if (!(row_array[start]==zero)){
1104        lastReducibleIndices[r]=start;
1105        return;
1106      }
1107    }
1108    lastReducibleIndices[r]=-1;
1109  }
1110  void backwardSubstitute(int r){
1111    int start=startIndices[r];
1112    assume(start<ncols);
1113    number_type zero=0;//npInit(0);
1114    number_type* row_array=rows[r];
1115    assume((!(npIsZero((number) row_array[start]))));
1116    assume(start<ncols);
1117    int other_row;
1118    if (!(npIsOne((number) row_array[r]))){
1119      //it should be one, but this safety is not expensive
1120      multiplyRow(r, F4mat_to_number_type(npInvers((number) row_array[start])));
1121    }
1122    int lastIndex=modP_lastIndexRow(row_array, ncols);
1123    assume(lastIndex<ncols);
1124    assume(lastIndex>=0);
1125    for(other_row=r-1;other_row>=0;other_row--){
1126      assume(lastReducibleIndices[other_row]<=start);
1127      if (lastReducibleIndices[other_row]==start){
1128        number_type* other_row_array=rows[other_row];
1129        number coef=npNeg((number) other_row_array[start]);
1130        assume(!(npIsZero(coef)));
1131        int i;
1132        assume(start>startIndices[other_row]);
1133        for(i=start;i<=lastIndex;i++){
1134          if (row_array[i]!=zero)
1135            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef,(number)row_array[i]),(number)other_row_array[i]));
1136        }
1137        updateLastReducibleIndex(other_row,r);
1138      }
1139    }
1140  }
1141  ~ModPMatrixBackSubstProxyOnArray<number_type>(){
1142    omfree(lastReducibleIndices);
1143  }
1144  void backwardSubstitute(){
1145    int i;
1146    for(i=nonZeroUntil;i>0;i--){
1147      backwardSubstitute(i);
1148    }
1149  }
1150};
1151template <class number_type > void simplest_gauss_modp(number_type* a, int nrows,int ncols){
1152  //use memmoves for changing rows
1153  if (TEST_OPT_PROT)
1154    PrintS("StartGauss\n");
1155  ModPMatrixProxyOnArray<number_type> mat(a,nrows,ncols);
1156 
1157  int c=0;
1158  int r=0;
1159  while(mat.findPivot(r,c)){
1160    //int pivot=find_pivot()
1161      mat.reduceOtherRowsForward(r);
1162    r++;
1163    c++;
1164  }
1165  ModPMatrixBackSubstProxyOnArray<number_type> backmat(mat);
1166  backmat.backwardSubstitute();
1167  //backward substitutions
1168  if (TEST_OPT_PROT)
1169    PrintS("StopGauss\n");
1170}
1171//int term_nodes_sort_crit(const void* a, const void* b);
1172template <class number_type> void noro_step(poly*p,int &pn,slimgb_alg* c){
1173  //Print("Input rows %d\n",pn);
1174  int j;
1175  if (TEST_OPT_PROT){
1176    Print("Input rows %d\n",pn);
1177  }
1178
1179  NoroCache<number_type> cache;
1180
1181  SparseRow<number_type> ** srows=(SparseRow<number_type>**) omalloc(pn*sizeof(SparseRow<number_type>*));
1182  int non_zeros=0;
1183  for(j=0;j<pn;j++){
1184   
1185    poly h=p[j];
1186    int h_len=pLength(h);
1187
1188    //number coef;
1189
1190
1191    srows[non_zeros]=noro_red_to_non_poly_t<number_type>(h,h_len,&cache,c);
1192    if (srows[non_zeros]!=NULL) non_zeros++;
1193  }
1194  std::vector<DataNoroCacheNode<number_type>*> irr_nodes;
1195  cache.collectIrreducibleMonomials(irr_nodes);
1196  //now can build up terms array
1197  //Print("historic irred Mon%d\n",cache.nIrreducibleMonomials);
1198  int n=irr_nodes.size();//cache.countIrreducibleMonomials();
1199  cache.nIrreducibleMonomials=n;
1200  if (TEST_OPT_PROT){
1201    Print("Irred Mon:%d\n",n);
1202    Print("red Mon:%d\n",cache.nReducibleMonomials);
1203  }
1204  TermNoroDataNode<number_type>* term_nodes=(TermNoroDataNode<number_type>*) omalloc(n*sizeof(TermNoroDataNode<number_type>));
1205 
1206  for(j=0;j<n;j++){
1207    assume(irr_nodes[j]!=NULL);
1208    assume(irr_nodes[j]->value_len==NoroCache<number_type>::backLinkCode);
1209    term_nodes[j].t=irr_nodes[j]->value_poly;
1210    assume(term_nodes[j].t!=NULL);
1211    term_nodes[j].node=irr_nodes[j];
1212  }
1213 
1214 
1215  qsort(term_nodes,n,sizeof(TermNoroDataNode<number_type>),term_nodes_sort_crit<number_type>);
1216  poly* terms=(poly*) omalloc(n*sizeof(poly));
1217
1218  int* old_to_new_indices=(int*) omalloc(cache.nIrreducibleMonomials*sizeof(int));
1219  for(j=0;j<n;j++){
1220    old_to_new_indices[term_nodes[j].node->term_index]=j;
1221    term_nodes[j].node->term_index=j;
1222    terms[j]=term_nodes[j].t;
1223  }
1224
1225  //if (TEST_OPT_PROT)
1226  //  Print("Evaluate Rows \n");
1227  pn=non_zeros;
1228  number_type* number_array=(number_type*) omalloc(n*pn*sizeof(number_type));
1229  memset(number_array,0,sizeof(number_type)*n*pn);
1230  number zero=npInit(0);
1231
1232  for(j=0;j<pn;j++){
1233    int i;
1234    number_type* row=number_array+n*j;
1235    /*for(i=0;i<n;i++){
1236      row[i]=zero;
1237    }*/
1238
1239    SparseRow<number_type>* srow=srows[j];
1240    if (srow){
1241    for(i=0;i<srow->len;i++){
1242      int idx=old_to_new_indices[srow->idx_array[i]];
1243      row[idx]=F4mat_to_number_type(srow->coef_array[i]);
1244    }
1245    delete srow;
1246    }
1247  }
1248 
1249  static int export_n=0;
1250  //export_mat(number_array,pn,n,"mat%i.py",++export_n);
1251  simplest_gauss_modp(number_array,pn,n);
1252
1253  int p_pos=0;
1254  for(j=0;j<pn;j++){
1255    poly h=row_to_poly(number_array+j*n,terms,n,c->r);
1256    if(h!=NULL){
1257      p[p_pos++]=h;
1258    }
1259  }
1260  pn=p_pos;
1261  omfree(terms);
1262  omfree(term_nodes);
1263  omfree(number_array);
1264  #ifdef NORO_NON_POLY
1265  omfree(srows);
1266  omfree(old_to_new_indices);
1267  #endif
1268  //don't forget the rank
1269 
1270}
1271
1272template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials( std::vector<DataNoroCacheNode<number_type> *>& res){
1273  int i;
1274  for(i=0;i<root.branches_len;i++){
1275    collectIrreducibleMonomials(1,root.branches[i],res);
1276  }
1277}
1278template <class number_type> void NoroCache<number_type>::collectIrreducibleMonomials(int level, NoroCacheNode* node, std::vector<DataNoroCacheNode<number_type>*>& res){
1279  assume(level>=0);
1280  if (node==NULL) return;
1281  if (level<pVariables){
1282    int i,sum;
1283    for(i=0;i<node->branches_len;i++){
1284      collectIrreducibleMonomials(level+1,node->branches[i],res);
1285    }
1286  } else {
1287    DataNoroCacheNode<number_type>* dn=(DataNoroCacheNode<number_type>*) node;
1288    if (dn->value_len==backLinkCode){
1289      res.push_back(dn);
1290    } 
1291  }
1292}
1293
1294template<class number_type> DataNoroCacheNode<number_type>* NoroCache<number_type>::getCacheReference(poly term){
1295  int i;
1296  NoroCacheNode* parent=&root;
1297  for(i=1;i<pVariables;i++){
1298    parent=parent->getBranch(p_GetExp(term,i,currRing));
1299    if (!(parent)){
1300      return NULL;
1301    }
1302  }
1303  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1304  return res_holder;
1305}
1306template<class number_type> poly NoroCache<number_type>::lookup(poly term, BOOLEAN& succ, int & len){
1307  int i;
1308  NoroCacheNode* parent=&root;
1309  for(i=1;i<pVariables;i++){
1310    parent=parent->getBranch(p_GetExp(term,i,currRing));
1311    if (!(parent)){
1312      succ=FALSE;
1313      return NULL;
1314    }
1315  }
1316  DataNoroCacheNode<number_type>* res_holder=(DataNoroCacheNode<number_type>*) parent->getBranch(p_GetExp(term,i,currRing));
1317  if (res_holder){
1318    succ=TRUE;
1319    if ((res_holder->value_len==backLinkCode)){
1320      len=1;
1321      return term;
1322    }
1323    len=res_holder->value_len;
1324    return res_holder->value_poly;
1325  } else {
1326    succ=FALSE;
1327    return NULL;
1328  }
1329}
1330#endif
1331
1332#endif
Note: See TracBrowser for help on using the repository browser.