source: git/kernel/tgb_internal.h @ 7b2896

jengelh-datetimespielwiese
Last change on this file since 7b2896 was 7b2896, checked in by Michael Brickenstein <bricken@…>, 16 years ago
*bricken: tuned noro F4 git-svn-id: file:///usr/local/Singular/svn/trunk@9887 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.1 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.58 2007-02-23 08:34:40 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};
99class DataNoroCacheNode;
100class 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};
119class MonRedResNP{
120public:
121  number coef;
122
123
124  DataNoroCacheNode* 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
140
141class NoroPlaceHolder{
142public:
143  DataNoroCacheNode* ref;
144  number coef;
145};
146
147//static ideal debug_Ideal;
148
149
150struct poly_list_node{
151  poly p;
152  poly_list_node* next;
153};
154
155struct int_pair_node{
156  int_pair_node* next;
157  int a;
158  int b;
159};
160struct monom_poly{
161  poly m;
162  poly f;
163};
164struct mp_array_list{
165  monom_poly* mp;
166  int size;
167  mp_array_list* next;
168};
169
170
171struct poly_array_list{
172  poly* p;
173  int size;
174  poly_array_list* next;
175};
176class slimgb_alg
177{
178  public:
179    slimgb_alg(ideal I, int syz_comp,BOOLEAN F4);
180                void introduceDelayedPairs(poly* pa,int s);
181    virtual ~slimgb_alg();
182    void cleanDegs(int lower, int upper);
183#ifndef HAVE_BOOST
184#ifdef USE_STDVECBOOL
185  vector<vector<bool> > states;
186#else
187  char** states;
188#endif
189#else
190  vector<dynamic_bitset<> > states;
191#endif
192  ideal add_later;
193  ideal S;
194  ring r;
195  int* lengths;
196  wlen_type* weighted_lengths;
197  long* short_Exps;
198  kStrategy strat;
199  int* T_deg;
200  int* T_deg_full;
201  poly tmp_lm;
202  poly* tmp_pair_lm;
203  sorted_pair_node** tmp_spn;
204  poly* expandS;
205  poly* gcd_of_terms;
206  int_pair_node* soon_free;
207  sorted_pair_node** apairs;
208  #if 0
209  BOOLEAN* modifiedS;
210  #endif
211  #ifdef TGB_RESORT_PAIRS
212  bool* replaced;
213  #endif
214  poly_list_node* to_destroy;
215  //for F4
216  mp_array_list* F;
217  poly_array_list* F_minus;
218
219  //end for F4
220#ifdef HEAD_BIN
221  struct omBin_s*   HeadBin;
222#endif
223  unsigned int reduction_steps;
224  int n;
225  //! array_lengths should be greater equal n;
226  int syz_comp;
227  int array_lengths; 
228  int normal_forms;
229  int current_degree;
230  int Rcounter;
231  int last_index;
232  int max_pairs;
233  int pair_top;
234  int easy_product_crit;
235  int extended_product_crit;
236  int average_length;
237  int lastDpBlockStart;
238  int lastCleanedDeg;
239  BOOLEAN isDifficultField;
240  BOOLEAN completed;
241  BOOLEAN is_homog;
242  BOOLEAN tailReductions;
243  BOOLEAN eliminationProblem;
244  BOOLEAN F4_mode;
245  BOOLEAN nc;
246  #ifdef TGB_RESORT_PAIRS
247  BOOLEAN used_b;
248  #endif
249};
250class red_object{
251 public:
252  kBucket_pt bucket;
253  poly p;
254  unsigned long sev;
255  int sugar;
256  void flatten();
257  void validate();
258  wlen_type initial_quality;
259  void adjust_coefs(number c_r, number c_ac_r);
260  wlen_type guess_quality(slimgb_alg* c);
261  int clear_to_poly();
262  void canonicalize();
263};
264
265
266enum calc_state
267  {
268    UNCALCULATED,
269    HASTREP//,
270    //UNIMPORTANT,
271    //SOONTREP
272  };
273static BOOLEAN pair_cmp(sorted_pair_node* a,sorted_pair_node* b);
274template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set);
275static int add_to_reductors(slimgb_alg* c, poly h, int len, int ecart, BOOLEAN simplified=FALSE);
276static int bucket_guess(kBucket* bucket);
277static poly redNFTail (poly h,const int sl,kStrategy strat, int len);
278static poly redNF2 (poly h,slimgb_alg* c , int &len, number&  m,int n=0);
279void free_sorted_pair_node(sorted_pair_node* s, ring r);
280static void shorten_tails(slimgb_alg* c, poly monom);
281static void replace_pair(int & i, int & j, slimgb_alg* c);
282//static sorted_pair_node** add_to_basis(poly h, int i, int j,slimgb_alg* c, int* ip=NULL);
283static void do_this_spoly_stuff(int i,int j,slimgb_alg* c);
284//ideal t_rep_gb(ring r,ideal arg_I);
285static BOOLEAN has_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* state);
286static int* make_connections(int from, poly bound, slimgb_alg* c);
287static int* make_connections(int from, int to, poly bound, slimgb_alg* c);
288void now_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
289static void soon_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c);
290static int pLcmDeg(poly a, poly b);
291static int simple_posInS (kStrategy strat, poly p,int len, wlen_type wlen);
292static BOOLEAN find_next_pair(slimgb_alg* c, BOOLEAN go_higher=TRUE);
293
294static sorted_pair_node* pop_pair(slimgb_alg* c);
295static BOOLEAN no_pairs(slimgb_alg* c);
296void clean_top_of_pair_list(slimgb_alg* c);
297static void super_clean_top_of_pair_list(slimgb_alg* c);
298static BOOLEAN state_is(calc_state state, const int & i, const int & j, slimgb_alg* c);
299static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, slimgb_alg* c=NULL);
300static int tgb_pair_better_gen(const void* ap,const void* bp);
301static poly redTailShort(poly h, kStrategy strat);
302static poly gcd_of_terms(poly p, ring r);
303static BOOLEAN extended_product_criterion(poly p1, poly gcd1, poly p2, poly gcd2, slimgb_alg* c);
304static poly kBucketGcd(kBucket* b, ring r);
305static void multi_reduction(red_object* los, int & losl, slimgb_alg* c);
306int slim_nsize(number n, ring r);
307sorted_pair_node* quick_pop_pair(slimgb_alg* c);
308sorted_pair_node* top_pair(slimgb_alg* c);
309sorted_pair_node** add_to_basis_ideal_quotient(poly h, slimgb_alg* c, int* ip);//, BOOLEAN new_pairs=TRUE);
310sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,slimgb_alg* c);
311int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj);
312int tgb_pair_better_gen2(const void* ap,const void* bp);
313int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev);
314//static int quality(poly p, int len, slimgb_alg* c);
315/**
316   makes on each red_object in a region a single_step
317 **/
318class reduction_step{
319 public:
320  /// we assume hat all occuring red_objects have same lm, and all
321  /// occ. lm's in r[l...u] are the same, only reductor does not occur
322  virtual void reduce(red_object* r, int l, int u);
323  //int reduction_id;
324  virtual ~reduction_step();
325  slimgb_alg* c;
326  int reduction_id;
327};
328class simple_reducer:public reduction_step{
329 public:
330  poly p;
331  kBucket_pt fill_back;
332  int p_len;
333  int reducer_deg;
334  simple_reducer(poly p, int p_len,int reducer_deg, slimgb_alg* c =NULL){
335    this->p=p;
336    this->reducer_deg=reducer_deg;
337    assume(p_len==pLength(p));
338    this->p_len=p_len;
339    this->c=c;
340  }
341  virtual void pre_reduce(red_object* r, int l, int u);
342  virtual void reduce(red_object* r, int l, int u);
343  ~simple_reducer();
344
345
346  virtual void do_reduce(red_object & ro);
347};
348
349//class sum_canceling_reducer:public reduction_step {
350//  void reduce(red_object* r, int l, int u);
351//};
352struct find_erg{
353  poly expand;
354  int expand_length;
355  int to_reduce_u;
356  int to_reduce_l;
357  int reduce_by;//index of reductor
358  BOOLEAN fromS;//else from los
359
360};
361
362static void multi_reduce_step(find_erg & erg, red_object* r, slimgb_alg* c);
363static void finalize_reduction_step(reduction_step* r);
364
365template <class len_type, class set_type>  int pos_helper(kStrategy strat, poly p, len_type len, set_type setL, polyset set){
366  //Print("POSHELER:%d",sizeof(wlen_type));
367  int length=strat->sl;
368  int i;
369  int an = 0;
370  int en= length;
371
372  if ((len>setL[length])
373      || ((len==setL[length]) && (pLmCmp(set[length],p)== -1)))
374    return length+1;
375
376  loop
377  {
378    if (an >= en-1)
379    {
380      if ((len<setL[an])
381          || ((len==setL[an]) && (pLmCmp(set[an],p) == 1))) return an;
382      return en;
383    }
384    i=(an+en) / 2;
385    if ((len<setL[i])
386        || ((len==setL[i]) && (pLmCmp(set[i],p) == 1))) en=i;
387    //else if ((len>setL[i])
388    //|| ((len==setL[i]) && (pLmCmp(set[i],p) == -1))) an=i;
389    else an=i;
390  }
391
392}
393#ifdef NORO_CACHE
394typedef unsigned short tgb_uint16;
395typedef unsigned char tgb_uint8;
396typedef unsigned int tgb_uint32;
397class NoroCacheNode{
398public:
399  NoroCacheNode** branches;
400  int branches_len;
401
402 
403  NoroCacheNode(){
404    branches=NULL;
405    branches_len=0;
406   
407  }
408  NoroCacheNode* setNode(int branch, NoroCacheNode* node){
409    if (branch>=branches_len){
410      if (branches==NULL){
411        branches_len=branch+1;
412        branches_len=si_max(branches_len,3);
413        branches=(NoroCacheNode**) omalloc(branches_len*sizeof(NoroCacheNode*));
414        int i;
415        for(i=0;i<branches_len;i++){
416          branches[i]=NULL;
417        }
418      }else{
419        int branches_len_old=branches_len;
420        branches_len=branch+1;
421        branches=(NoroCacheNode**) omrealloc(branches,branches_len*sizeof(NoroCacheNode*));
422        int i;
423        for(i=branches_len_old;i<branches_len;i++){
424          branches[i]=NULL;
425        }
426      }
427    }
428    assume(branches[branch]==NULL);
429    branches[branch]=node;
430    return node;
431  }
432  NoroCacheNode* getBranch(int branch){
433    if (branch<branches_len) return branches[branch];
434    return NULL;
435  }
436  virtual ~NoroCacheNode(){
437    int i;
438    for(i=0;i<branches_len;i++){
439      delete branches[i];
440    }
441    omfree(branches);
442  }
443  NoroCacheNode* getOrInsertBranch(int branch){
444    if ((branch<branches_len)&&(branches[branch]))
445      return branches[branch];
446    else{
447      return setNode(branch,new NoroCacheNode());
448    }
449  }
450};
451class DenseRow{
452public:
453  number* array;
454  int begin;
455  int end;
456  DenseRow(){
457    array=NULL;
458  }
459  ~DenseRow(){
460    omfree(array);
461  }
462};
463class SparseRow{
464public:
465  int* idx_array;
466  number* coef_array;
467  int len;
468  SparseRow(){
469    len=0;
470    idx_array=NULL;
471    coef_array=NULL;
472  }
473  SparseRow(int n){
474    len=n;
475    idx_array=(int*) omalloc(n*sizeof(int));
476    coef_array=(number*) omalloc(n*sizeof(number));
477  }
478  ~SparseRow(){
479    omfree(idx_array);
480    omfree(coef_array);
481  }
482};
483
484class DataNoroCacheNode:public NoroCacheNode{
485public:
486 
487  int value_len;
488  poly value_poly;
489  #ifdef NORO_SPARSE_ROWS_PRE
490  SparseRow* row;
491  #else
492  DenseRow* row;
493  #endif
494  int term_index;
495  DataNoroCacheNode(poly p, int len){
496    value_len=len;
497    value_poly=p;
498    row=NULL;
499    term_index=-1;
500  }
501  #ifdef NORO_SPARSE_ROWS_PRE
502  DataNoroCacheNode(SparseRow* row){
503    if (row!=NULL)
504      value_len=row->len;
505    else
506      value_len=0;
507    value_poly=NULL;
508    this->row=row;
509    term_index=-1;
510  }
511  #endif
512  ~DataNoroCacheNode(
513  ){
514    //p_Delete(&value_poly,currRing);
515    if (row) delete row;
516  }
517};
518class NoroCache{
519public:
520  poly temp_term;
521  void evaluatePlaceHolder(number* row,std::vector<NoroPlaceHolder>& place_holders);
522  void collectIrreducibleMonomials( std::vector<DataNoroCacheNode*>& res);
523  void collectIrreducibleMonomials(int level,  NoroCacheNode* node, std::vector<DataNoroCacheNode*>& res);
524  void evaluateRows();
525  void evaluateRows(int level, NoroCacheNode* node);
526#ifdef NORO_RED_ARRAY_RESERVER
527  int reserved;
528  poly* recursionPolyBuffer;
529#endif
530  static const int backLinkCode=-222;
531  DataNoroCacheNode* insert(poly term, poly nf, int len){
532    //assume(impl.find(p_Copy(term,currRing))==impl.end());
533    //assume(len==pLength(nf));
534    assume(npIsOne(p_GetCoeff(term,currRing)));
535    if (term==nf){
536      term=p_Copy(term,currRing);
537
538      ressources.push_back(term);
539      nIrreducibleMonomials++;
540      return treeInsertBackLink(term);
541     
542    } else{
543     
544      if (nf){
545        //nf=p_Copy(nf,currRing);
546        assume(p_LmCmp(nf,term,currRing)==-1);
547        ressources.push_back(nf);
548      }
549      return treeInsert(term,nf,len);
550     
551    }
552   
553    //impl[term]=std::pair<PolySimple,int> (nf,len);
554  }
555  #ifdef NORO_SPARSE_ROWS_PRE
556  DataNoroCacheNode* insert(poly term, SparseRow* srow){
557    //assume(impl.find(p_Copy(term,currRing))==impl.end());
558    //assume(len==pLength(nf));
559
560      return treeInsert(term,srow);
561     
562 
563    //impl[term]=std::pair<PolySimple,int> (nf,len);
564  }
565  #endif
566  DataNoroCacheNode* insertAndTransferOwnerShip(poly t, ring r){
567   
568    ressources.push_back(t);
569    DataNoroCacheNode* res=treeInsertBackLink(t);
570    res->term_index=nIrreducibleMonomials;
571    nIrreducibleMonomials++;
572    return res;
573  }
574  poly lookup(poly term, BOOLEAN& succ, int & len);
575  DataNoroCacheNode* getCacheReference(poly term);
576  NoroCache(){
577    buffer=NULL;
578#ifdef NORO_RED_ARRAY_RESERVER
579    reserved=0;
580    recursionPolyBuffer=(poly*)omalloc(1000000*sizeof(poly));
581#endif
582    nIrreducibleMonomials=0;
583    nReducibleMonomials=0;
584    temp_term=pOne();
585    tempBufferSize=3000;
586    tempBuffer=omalloc(tempBufferSize);
587  }
588  void ensureTempBufferSize(size_t size){
589    if (tempBufferSize<size){
590      tempBufferSize=2*size;
591      omfree(tempBuffer);
592      tempBuffer=omalloc(tempBufferSize);
593    }
594  }
595#ifdef NORO_RED_ARRAY_RESERVER
596  poly* reserve(int n){
597    poly* res=recursionPolyBuffer+reserved;
598    reserved+=n;
599    return res;
600  }
601  void free(int n){
602    reserved-=n;
603  }
604#endif
605  ~NoroCache(){
606    int s=ressources.size();
607    int i;
608    for(i=0;i<s;i++){
609      p_Delete(&ressources[i].impl,currRing);
610    }
611    p_Delete(&temp_term,currRing);
612#ifdef NORO_RED_ARRAY_RESERVER
613    omfree(recursionPolyBuffer);
614#endif
615   omfree(tempBuffer);
616  }
617 
618  int nIrreducibleMonomials;
619  int nReducibleMonomials;
620  void* tempBuffer;
621  size_t tempBufferSize;
622protected:
623  DataNoroCacheNode* treeInsert(poly term,poly nf,int len){
624    int i;
625    nReducibleMonomials++;
626    int nvars=pVariables;
627    NoroCacheNode* parent=&root;
628    for(i=1;i<nvars;i++){
629      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
630    }
631    return (DataNoroCacheNode*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode(nf,len));
632  }
633  #ifdef NORO_SPARSE_ROWS_PRE
634  DataNoroCacheNode* treeInsert(poly term,SparseRow* srow){
635    int i;
636    nReducibleMonomials++;
637    int nvars=pVariables;
638    NoroCacheNode* parent=&root;
639    for(i=1;i<nvars;i++){
640      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
641    }
642    return (DataNoroCacheNode*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode(srow));
643  }
644  #endif
645  DataNoroCacheNode* treeInsertBackLink(poly term){
646    int i;
647    int nvars=pVariables;
648    NoroCacheNode* parent=&root;
649    for(i=1;i<nvars;i++){
650      parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing));
651    }
652    return (DataNoroCacheNode*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode(term,backLinkCode));
653  }
654 
655  //@TODO descruct nodes;
656  typedef std::vector<PolySimple> poly_vec;
657  poly_vec ressources;
658  //typedef std::map<PolySimple,std::pair<PolySimple,int> > cache_map;
659  //cache_map impl;
660  NoroCacheNode root;
661  number* buffer;
662};
663MonRedResNP noro_red_mon_to_non_poly(poly t,  NoroCache* cache,slimgb_alg* c);
664template<class storage_type> SparseRow* noro_red_to_non_poly_t(poly p, int &len, NoroCache* cache,slimgb_alg* c){
665  assume(len==pLength(p));
666  poly orig_p=p;
667  if (p==NULL) {
668    len=0;
669    return NULL;
670  }
671 
672  number zero=npInit(0);
673  MonRedResNP* mon=(MonRedResNP*) omalloc(len*sizeof(MonRedResNP));
674  int i=0;
675
676  while(p){
677
678    poly t=p;
679    pIter(p);
680    pNext(t)=NULL;
681   
682#ifndef NDEBUG
683    number coef_debug=p_GetCoeff(t,currRing);
684#endif
685    MonRedResNP red=noro_red_mon_to_non_poly(t,cache,c);
686    mon[i]=red;
687    i++;
688  }
689 
690  assume(i==len);
691  len=i;
692  //in the loop before nIrreducibleMonomials increases, so position here is important
693  size_t temp_size_bytes=cache->nIrreducibleMonomials*sizeof(storage_type);
694  cache->ensureTempBufferSize(temp_size_bytes);
695  storage_type* temp_array=(storage_type*) cache->tempBuffer;//omalloc(cache->nIrreducibleMonomials*sizeof(storage_type));
696  int temp_size=cache->nIrreducibleMonomials;
697  memset(temp_array,0,temp_size_bytes);
698  for(i=0;i<len;i++){
699    MonRedResNP red=mon[i];
700    if ((red.ref)){
701      if (red.ref->row){
702        SparseRow* row=red.ref->row;
703        number coef=red.coef;
704        int j;
705        if (!(npIsOne(coef))){
706        for(j=0;j<row->len;j++){
707          int idx=row->idx_array[j];
708          assume(!(npIsZero(coef)));
709          assume(!(npIsZero(row->coef_array[j])));
710          temp_array[idx]=(storage_type) (unsigned int) npAddM((number) temp_array[idx],npMultM(row->coef_array[j],coef));
711          assume(idx<temp_size);
712        }}else{
713          for(j=0;j<row->len;j++){
714            int idx=row->idx_array[j];
715            temp_array[idx]=(storage_type) (unsigned int) npAddM((number) temp_array[idx],row->coef_array[j]);
716            assume(idx<temp_size);
717          }
718        }
719      }
720      else{
721        if (red.ref->value_len==NoroCache::backLinkCode){
722          temp_array[red.ref->term_index]=(storage_type) (unsigned int) npAddM((number) temp_array[red.ref->term_index],red.coef);
723        } else {
724          //PrintS("third case\n");
725        }
726      }
727    }
728  }
729  int non_zeros=0;
730  for(i=0;i<cache->nIrreducibleMonomials;i++){
731    //if (!(temp_array[i]==0)){
732    //  non_zeros++;
733    //}
734    assume(((temp_array[i]!=0)==0)|| (((temp_array[i]!=0)==1)));
735    non_zeros+=(temp_array[i]!=0);
736  }
737 
738  if (non_zeros==0){
739    omfree(mon);
740    return NULL;
741  }
742  SparseRow* res=new SparseRow(non_zeros);
743  int pos=0;
744  for(i=0;i<cache->nIrreducibleMonomials;i++){
745    if (!(0==temp_array[i])){
746   
747      res->idx_array[pos]=i;
748      res->coef_array[pos]=(number) temp_array[i];
749
750      pos++;
751      non_zeros--;
752      if (non_zeros==0) break;
753    }
754   
755  }
756  //omfree(temp_array);
757
758  omfree(mon);
759  return res;
760}
761#endif
762static wlen_type pair_weighted_length(int i, int j, slimgb_alg* c);
763wlen_type pELength(poly p, ring r);
764int terms_sort_crit(const void* a, const void* b);
765//void simplest_gauss_modp(number* a, int nrows,int ncols);
766// a: a[0,0],a[0,1]....a[nrows-1,ncols-1]
767// assume: field is Zp
768#ifdef USE_NORO
769#define slim_prec_cast(a) (unsigned int) (a)
770#define F4mat_to_number_type(a) (number_type) slim_prec_cast(a)
771
772template <class number_type > void write_poly_to_row(number_type* row, poly h, poly*terms, int tn, ring r){
773  //poly* base=row;
774  while(h!=NULL){
775    //Print("h:%i\n",h);
776    number coef=p_GetCoeff(h,r);
777    poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit);
778    assume(ptr_to_h!=NULL);
779    int pos=ptr_to_h-terms;
780    row[pos]=F4mat_to_number_type(coef);
781    //number_type_array[base+pos]=coef;
782    pIter(h);
783  }
784}
785template <class number_type > poly row_to_poly(number_type* row, poly* terms, int tn, ring r){
786  poly h=NULL;
787  int j;
788  number_type zero=0;//;npInit(0);
789  for(j=tn-1;j>=0;j--){
790    if (!(zero==(row[j]))){
791      poly t=terms[j];
792      t=p_LmInit(t,r);
793      p_SetCoeff(t,(number) row[j],r);
794      pNext(t)=h;
795      h=t;
796    }
797   
798  }
799  return h;
800}
801template <class number_type > int modP_lastIndexRow(number_type* row,int ncols){
802  int lastIndex;
803  const number_type zero=0;//npInit(0);
804  for(lastIndex=ncols-1;lastIndex>=0;lastIndex--){
805    if (!(row[lastIndex]==zero)){
806      return lastIndex;
807    }
808  }
809  return -1;
810}
811template <class number_type>class ModPMatrixBackSubstProxyOnArray;
812template <class number_type > class ModPMatrixProxyOnArray{
813public:
814  friend class ModPMatrixBackSubstProxyOnArray<number_type>;
815               
816  int ncols,nrows;
817  ModPMatrixProxyOnArray(number_type* array, int nrows, int ncols){
818    this->ncols=ncols;
819    this->nrows=nrows;
820    rows=(number_type**) omalloc(nrows*sizeof(number_type*));
821    startIndices=(int*)omalloc(nrows*sizeof(int));
822    int i;
823    for(i=0;i<nrows;i++){
824      rows[i]=array+(i*ncols);
825      updateStartIndex(i,-1);
826    }
827  }
828  ~ModPMatrixProxyOnArray(){
829    omfree(rows);
830    omfree(startIndices);
831  }
832 
833  void permRows(int i, int j){
834    number_type* h=rows[i];
835    rows[i]=rows[j];
836    rows[j]=h;
837    int hs=startIndices[i];
838    startIndices[i]=startIndices[j];
839    startIndices[j]=hs;
840  }
841  void multiplyRow(int row, number_type coef){
842    int i;
843    number_type* row_array=rows[row];
844    for(i=startIndices[row];i<ncols;i++){
845      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
846    }
847  }
848  void reduceOtherRowsForward(int r){
849
850    //assume rows "under r" have bigger or equal start index
851    number_type* row_array=rows[r];
852    number_type zero=F4mat_to_number_type(npInit(0));
853    int start=startIndices[r];
854    number_type coef=row_array[start];
855    assume(start<ncols);
856    int other_row;
857    assume(!(npIsZero((number) row_array[start])));
858    if (!(npIsOne((number) coef)))
859      multiplyRow(r,F4mat_to_number_type(npInvers((number) coef)));
860    assume(npIsOne((number) row_array[start]));
861    int lastIndex=modP_lastIndexRow(row_array, ncols);
862    number minus_one=npInit(-1);
863    for (other_row=r+1;other_row<nrows;other_row++){
864      assume(startIndices[other_row]>=start);
865      if (startIndices[other_row]==start){
866        int i;
867        number_type* other_row_array=rows[other_row];
868        number coef2=npNeg((number) other_row_array[start]);
869        if (coef2==minus_one){
870          for(i=start;i<=lastIndex;i++){
871            if (row_array[i]!=zero)
872              other_row_array[i]=F4mat_to_number_type(npSubM((number) other_row_array[i], (number) row_array[i]));
873          }
874      }else {
875          //assume(FALSE);
876          for(i=start;i<=lastIndex;i++){
877            if (row_array[i]!=zero)
878            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef2,(number) row_array[i]),(number) other_row_array[i]));
879          }
880        }
881        updateStartIndex(other_row,start);
882        assume(npIsZero((number) other_row_array[start]));
883      }
884    }
885  }
886  void updateStartIndex(int row,int lower_bound){
887    number_type* row_array=rows[row];
888    assume((lower_bound<0)||(npIsZero((number) row_array[lower_bound])));
889    int i;
890    //number_type zero=npInit(0);
891    for(i=lower_bound+1;i<ncols;i++){
892      if (!(row_array[i]==0))
893        break;
894    }
895    startIndices[row]=i;
896  }
897  int getStartIndex(int row){
898    return startIndices[row];
899  }
900  BOOLEAN findPivot(int &r, int &c){
901    //row>=r, col>=c
902   
903    while(c<ncols){
904      int i;
905      for(i=r;i<nrows;i++){
906        assume(startIndices[i]>=c);
907        if (startIndices[i]==c){
908          //r=i;
909          if (r!=i)
910            permRows(r,i);
911          return TRUE;
912        }
913      }
914      c++;
915    }
916    return FALSE;
917  }
918protected:
919  number_type** rows;
920  int* startIndices;
921};
922template <class number_type > class ModPMatrixBackSubstProxyOnArray{
923  int *startIndices;
924  number_type** rows;
925  int *lastReducibleIndices;
926  int ncols;
927  int nrows;
928  int nonZeroUntil;
929public:
930  void multiplyRow(int row, number_type coef){
931    int i;
932    number_type* row_array=rows[row];
933    for(i=startIndices[row];i<ncols;i++){
934      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
935    }
936  }
937  ModPMatrixBackSubstProxyOnArray<number_type> (ModPMatrixProxyOnArray<number_type> & p){
938//  (number_type* array, int nrows, int ncols, int* startIndices, number_type** rows){
939    //we borrow some parameters ;-)
940    //we assume, that nobody changes the order of the rows
941    this->startIndices=p.startIndices;
942    this->rows=p.rows;
943    this->ncols=p.ncols;
944    this->nrows=p.nrows;
945    lastReducibleIndices=(int*) omalloc(nrows*sizeof(int));
946    nonZeroUntil=0;
947    while(nonZeroUntil<nrows){
948      if (startIndices[nonZeroUntil]<ncols){
949       
950        nonZeroUntil++;
951      } else break;
952     
953    }
954    if (TEST_OPT_PROT)
955      Print("rank:%i\n",nonZeroUntil);
956    nonZeroUntil--;
957    int i;
958    for(i=0;i<=nonZeroUntil;i++){
959      assume(startIndices[i]<ncols);
960      assume(!(npIsZero((number) rows[i][startIndices[i]])));
961      assume(startIndices[i]>=i);
962      updateLastReducibleIndex(i,nonZeroUntil+1);
963    }
964  }
965  void updateLastReducibleIndex(int r, int upper_bound){
966    number_type* row_array=rows[r];
967    if (upper_bound>nonZeroUntil) upper_bound=nonZeroUntil+1;
968    int i;
969    const number_type zero=0;//npInit(0);
970    for(i=upper_bound-1;i>r;i--){
971      int start=startIndices[i];
972      assume(start<ncols);
973      if (!(row_array[start]==zero)){
974        lastReducibleIndices[r]=start;
975        return;
976      }
977    }
978    lastReducibleIndices[r]=-1;
979  }
980  void backwardSubstitute(int r){
981    int start=startIndices[r];
982    assume(start<ncols);
983    number_type zero=0;//npInit(0);
984    number_type* row_array=rows[r];
985    assume((!(npIsZero((number) row_array[start]))));
986    assume(start<ncols);
987    int other_row;
988    if (!(npIsOne((number) row_array[r]))){
989      //it should be one, but this safety is not expensive
990      multiplyRow(r, F4mat_to_number_type(npInvers((number) row_array[start])));
991    }
992    int lastIndex=modP_lastIndexRow(row_array, ncols);
993    assume(lastIndex<ncols);
994    assume(lastIndex>=0);
995    for(other_row=r-1;other_row>=0;other_row--){
996      assume(lastReducibleIndices[other_row]<=start);
997      if (lastReducibleIndices[other_row]==start){
998        number_type* other_row_array=rows[other_row];
999        number coef=npNeg((number) other_row_array[start]);
1000        assume(!(npIsZero(coef)));
1001        int i;
1002        assume(start>startIndices[other_row]);
1003        for(i=start;i<=lastIndex;i++){
1004          if (row_array[i]!=zero)
1005            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef,(number)row_array[i]),(number)other_row_array[i]));
1006        }
1007        updateLastReducibleIndex(other_row,r);
1008      }
1009    }
1010  }
1011  ~ModPMatrixBackSubstProxyOnArray<number_type>(){
1012    omfree(lastReducibleIndices);
1013  }
1014  void backwardSubstitute(){
1015    int i;
1016    for(i=nonZeroUntil;i>0;i--){
1017      backwardSubstitute(i);
1018    }
1019  }
1020};
1021template <class number_type > void simplest_gauss_modp(number_type* a, int nrows,int ncols){
1022  //use memmoves for changing rows
1023  if (TEST_OPT_PROT)
1024    PrintS("StartGauss\n");
1025  ModPMatrixProxyOnArray<number_type> mat(a,nrows,ncols);
1026 
1027  int c=0;
1028  int r=0;
1029  while(mat.findPivot(r,c)){
1030    //int pivot=find_pivot()
1031      mat.reduceOtherRowsForward(r);
1032    r++;
1033    c++;
1034  }
1035  ModPMatrixBackSubstProxyOnArray<number_type> backmat(mat);
1036  backmat.backwardSubstitute();
1037  //backward substitutions
1038  if (TEST_OPT_PROT)
1039    PrintS("StopGauss\n");
1040}
1041#endif
1042
1043#endif
Note: See TracBrowser for help on using the repository browser.