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:
|