File: gb-test1.hpp

package info (click to toggle)
macaulay2 1.21%2Bds-3
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 133,096 kB
  • sloc: cpp: 110,377; ansic: 16,306; javascript: 4,193; makefile: 3,821; sh: 3,580; lisp: 764; yacc: 590; xml: 177; python: 140; perl: 114; lex: 65; awk: 3
file content (397 lines) | stat: -rw-r--r-- 11,603 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
/* Copyright 2003, Michael E. Stillman */

#ifndef _gbB_h_
#define _gbB_h_

#include "comp-gb.hpp"

#include "gbring.hpp"
#include "montable.hpp"
#include "reducedgb.hpp"
#include "interreduce.hpp"
#include "f4/hilb-fcn.hpp"
#include "polyring.hpp"

// This GB algorithm does NOT handle:
//   coeffs == ZZ
//   local GB

class GBWeight;

const int ELEMB_IN_RING = 1;
const int ELEMB_MINGEN = 2;
const int ELEMB_MINGB = 4;

/**
    @ingroup gb

    @brief A modification of the default algorithm to try to speed it up in
   inhomogeneous cases.
*/
class gbB : public GBComputation
{
 public:
  typedef int gbelem_type;
  // ELEMB_IN_RING: a bit set.  If this is set, then the next two bits need to
  //   off.
  // ELEMB_MINGEN: a bit set: 0 means it is provably not needed to
  //   generated the ideal or submodule.  1 means it might be a minimal
  //   generator.
  //   NOTE: often the min gens obtained this way are not so good.  Exceptions
  //   are
  //    for graded submodules.  Then the meaning of 1 is: it IS a min gen.
  // ELEMB_MINGB: a bit set: 0 means is not part of the min gb (so far).  1
  // means it is.
  // SO: possible values are 0, 1, 2, 4, 6.
  // Test using code like this:
  //   if (g->minlevel & ELEMB_MINGEN) { .. do this if it is a possible min gen
  //   .. } else { .. if not .. };

  struct gbelem
  {
    POLY g;
    int deg;
    int gap;  // the homogenizing degree, called "alpha", or the "ecart".  It is
              // deg g - deg lead term of g.
    int reduced_deg;  // the poly might be of actual degree lower than 'deg'.
                      // This is the actual deg.
    int reduced_gap;  // This is the gap value assuming the degree is
                      // reduced_deg.  Used in finding divisors.
    int size;  // number of monomials, when the element is first inserted. After
    // auto reduction, the number can change.  It is not yet clear if we should
    // change this value...
    exponents lead;  // -1..nvars-1, the -1 part is the component
    gbelem_type minlevel;
  };

  /* Types of minimality */
  enum spair_type {
    SPAIR_SPAIR,
    SPAIR_RING,
    SPAIR_SKEW,
    SPAIR_GEN,
    SPAIR_ELEM
  };

 public:
  // This is only public to allow spair_sorter to use it!!
  struct spair
  {
    spair *next;
    spair_type type; /* SPAIR_SPAIR,
                        SPAIR_GEN, SPAIR_ELEM, SPAIR_RING, SPAIR_SKEW */
    int deg;
    exponents lcm; /* Contains homogenizing variable, component */
    union
    {
      POLY f; /* SPAIR_GEN, SPAIR_ELEM */
      struct pair
      {
        int i, j; /* i refers to a GB element. j refers to GB element
                     (SPAIR_SPAIR)
                     or a ring element (SPAIR_RING) or a variable number
                     (SPAIR_SKEW)
                   */
      } pair;
    } x;
    gbvector *lead_of_spoly;  // experimental
    gbvector *&f() { return x.f.f; }
    gbvector *&fsyz() { return x.f.fsyz; }
  };

 private:
  typedef VECTOR(spair *) spairs;

  struct SPairSet : public our_new_delete
  {
    int nelems;      /* Includes the number in this_set */
    int n_in_degree; /* The number in 'this_set */
    spair *heap;
    int n_computed; /* The number removed via next() */

    // Each of the following is initialized on each call to
    // spair_set_prepare_next_degree
    spair *spair_list;
    spair spair_deferred_list;   // list header
    spair *spair_last_deferred;  // points to last elem of previous list,
                                 // possibly the header

    spair *gen_list;
    spair gen_deferred_list;   // list header
    spair *gen_last_deferred;  // points to last elem of previous list, possibly
                               // the header

    SPairSet();
  };

 private:
  // Stashes
  stash *spair_stash;
  stash *gbelem_stash;
  size_t exp_size;  // in bytes
  stash *lcm_stash;

  // Data
  const PolynomialRing *originalR;
  GBRing *R;
  const GBWeight *weightInfo;
  M2_arrayint gb_weights;

  // Ring information

  const FreeModule *F;
  const FreeModule *Fsyz;
  int nvars;

  VECTOR(gbelem *) gb;  // Contains any quotient ring elements

  ReducedGB *minimal_gb;
  bool minimal_gb_valid;

  MonomialTable *ringtable;  // At most one of these two will be non-NULL.

  MonomialTable *lookup;

  SPairSet *S;

  Interreducer *G;

  VECTOR(gbvector *) _syz;

  int _strategy;
  int n_rows_per_syz;
  bool _collect_syz;
  bool _is_ideal;
  int first_gb_element; /* First index past the skew variable squares, quotient
                           elements,
                           in the array in G */
  int first_in_degree;  /* for the current 'sugar' degree, this is the first GB
                           element
                            inserted for this degree */
  int complete_thru_this_degree; /* Used in reporting status to the user */

  /* (new) state machine */
  enum gbB_state {
    STATE_NEWDEGREE,
    STATE_NEWPAIRS,
    STATE_SPAIRS,
    STATE_GENS,
    STATE_AUTOREDUCE,
    STATE_DONE
  } state;
  int this_degree;
  int npairs;  // in this_degree
  int np_i;
  int ar_i;
  int ar_j;
  int n_gb;

  /* stats */
  int stats_ngcd1;
  int stats_nreductions;
  int stats_ntail;
  int stats_ngb;
  int stats_npairs;

  /* for ending conditions */
  int n_syz;
  int n_pairs_computed;
  int n_gens_left;
  int n_subring;

  // Hilbert function information
  HilbertController *hilbert;  // null if Hilbert function not being used
  int n_saved_hilb;
#if 0
  bool use_hilb;
  bool hilb_new_elems;  // True if any new elements since HF was last computed
  int hilb_n_in_degree; // The number of new elements that we expect to find
                         // in this degree.
  const RingElement *hf_orig;   // The Hilbert function that we are given at the beginning
  RingElement *hf_diff;         // The difference between hf_orig and the computed hilb fcn
#endif
  //////////////////////////////////////////

  // Reduction count: used to defer spairs which are likely to reduce to 0
  long max_reduction_count;
  void spair_set_defer(spair *&p);

  // Stashing previous divisors;
  int divisor_previous;
  int divisor_previous_comp;

 private:
  exponents exponents_make();

  /* initialization */
  void initialize(const Matrix *m,
                  int csyz,
                  int nsyz,
                  M2_arrayint gb_weights,
                  int strat,
                  int max_reduction_count0);
  spair *new_gen(int i, gbvector *f, ring_elem denom);

  int gbelem_COMPONENT(gbelem *g) { return g->g.f->comp; }
  int spair_COMPONENT(spair *s)
  {
    // Only valid if this is an SPAIR_ELEM, SPAIR_RING, SPAIR_SKEW.
    // Probably better is to put it into spair structure.
    return gb[s->x.pair.i]->g.f->comp;
  }

  /* new state machine routines */
  void insert_gb(POLY f, gbelem_type minlevel);
  bool process_spair(spair *p);
  void do_computation();

  gbelem *gbelem_ring_make(gbvector *f);

  gbelem *gbelem_make(gbvector *f,     // grabs f
                      gbvector *fsyz,  // grabs fsyz
                      gbelem_type minlevel,
                      int deg);

  void gbelem_text_out(buffer &o, int i, int nterms = 3) const;

  /* spair creation */
  /* negative indices index quotient ring elements */
  spair *spair_node();
  spair *spair_make(int i, int j);
  spair *spair_make_gen(POLY f);
  spair *spair_make_skew(int i, int v);
  spair *spair_make_ring(int i, int j);
  void spair_text_out(buffer &o, spair *p);
  void spair_delete(spair *&p);

  /* spair handling */
  bool pair_not_needed(spair *p, gbelem *m);
  void remove_unneeded_pairs(int id);
  bool is_gcd_one_pair(spair *p);
  spairs::iterator choose_pair(spairs::iterator first, spairs::iterator next);
  void minimalize_pairs(spairs &new_set);
  void update_pairs(int id);

  /* spair set handling */
  void remove_spair_list(spair *&set);
  void remove_SPairSet();
  void spair_set_insert(spair *p);
  /* Insert a LIST of s pairs into S */
  spair *spair_set_next();
  /* Removes the next element of the current degree, returning NULL if none left
   */
  int spair_set_determine_next_degree(int &nextdegree);
  int spair_set_prepare_next_degree(int &nextdegree);
  /* Finds the next degree to consider, returning the number of spairs in that
   * degree */
  void spair_set_show_mem_usage();

  void spairs_sort(int len, spair *&list);
  void spairs_reverse(spair *&ps);

  void spair_set_lead_spoly(spair *p);

  /* Sorts the list of spairs 'list' (which has length 'len') */

  /* Hilbert function handling */
  void
  flush_pairs();  // Used to flush the rest of the pairs in the current degree.
  RingElement /* or null */ *
  compute_hilbert_function();  // Compute the HF of _hilb_matrix.
#if 0
  //TODO: remove once hilbert is functional
  Matrix *make_lead_term_matrix(); // The submodule of all lead terms
#endif

  /* reduction */
  void auto_reduce_by(int id);
  void compute_s_pair(spair *p);
  bool reduceit(spair *p);
  void collect_syzygy(gbvector *fsyz);

  enum ComputationStatusCode computation_is_complete();

  virtual bool stop_conditions_ok() { return true; }
 private:
  /* Making the minimal GB */

  void minimalize_gb();

  int find_good_divisor(exponents e, int x, int degf, int &result_gap);

  void remainder(POLY &f, int degf, bool use_denom, ring_elem &denom);
  // denom must start out as an element of the base R->getCoefficients().
  // denom is multiplied by the coefficient which is multiplied to f in order to
  // reduce f.
  // i.e. the result g satisfies: g = c*f mod GB, where new(denom) = c *
  // old(denom).

 public:
  //////////////////////////
  // Computation routines //
  //////////////////////////

  static gbB *create(const Matrix *m,
                     M2_bool collect_syz,
                     int n_rows_to_keep,
                     M2_arrayint gb_weights,
                     int strategy,
                     M2_bool use_max_degree,
                     int max_degree,
                     int max_reduction_count);

  static gbB *create_forced(const Matrix *m,
                            const Matrix *gb,
                            const Matrix *mchange);

  void remove_gb();
  virtual ~gbB();

  virtual int kind() { return 231; }  // FIX THIS!!
  void start_computation();

  virtual const PolynomialRing *get_ring() const { return originalR; }
  virtual Computation /* or null */ *set_hilbert_function(const RingElement *h);

  virtual const Matrix /* or null */ *get_gb();

  virtual const Matrix /* or null */ *get_mingens();

  virtual const Matrix /* or null */ *get_change();

  virtual const Matrix /* or null */ *get_syzygies();

  virtual const Matrix /* or null */ *get_initial(int nparts);

  virtual const Matrix /* or null */ *get_parallel_lead_terms(M2_arrayint w);

  virtual const Matrix /* or null */ *matrix_remainder(const Matrix *m);

  virtual M2_bool matrix_lift(const Matrix *m,
                              const Matrix /* or null */ **result_remainder,
                              const Matrix /* or null */ **result_quotient);

  virtual int contains(const Matrix *m);

  virtual void text_out(buffer &o) const;
  /* This displays statistical information, and depends on the
     M2_gbTrace value */

  virtual int complete_thru_degree() const;
  // The computation is complete up through this degree.

  /* Debug display routines */
  void debug_spair(spair *p);
  void debug_spairs(spair *spairlist);
  void debug_spair_array(spairs &spairlist);
  void show() const;

  void show_mem_usage();
};

#endif

// Local Variables:
// compile-command: "make -C $M2BUILDDIR/Macaulay2/e "
// indent-tabs-mode: nil
// End: