source: git/kernel/GBEngine/syz4.cc @ 0612fd

spielwiese
Last change on this file since 0612fd was 0612fd, checked in by Hans Schoenemann <hannes@…>, 5 years ago
unified (un)lukely/(UN)LIKELY -> (UN)LIKELY
  • Property mode set to 100644
File size: 24.1 KB
Line 
1/***************************************
2 *  Computer Algebra System SINGULAR   *
3 ***************************************/
4/*
5 * ABSTRACT: resolutions
6 * reference: https://arxiv.org/abs/1502.01654
7 */
8
9#include <kernel/GBEngine/syz.h>
10#include <omalloc/omalloc.h>
11#include <coeffs/numbers.h>
12#include <kernel/polys.h>
13#include <kernel/ideals.h>
14
15#include <vector>
16#include <map>
17
18/*
19 * set variables[i] to false if the i-th variable does not appear among the
20 * leading terms of L
21 */
22static void update_variables(std::vector<bool> &variables, const ideal L)
23{
24    const ring R = currRing;
25    const int l = L->ncols-1;
26    int k;
27    for (int j = R->N; j > 0; j--)
28    {
29        if (variables[j-1])
30        {
31            for (k = l; k >= 0; k--)
32            {
33                if (p_GetExp(L->m[k], j, R) > 0)
34                {
35                    break;
36                }
37            }
38            if (k < 0)
39            {   // no break
40                variables[j-1] = false;
41            }
42        }
43    }
44}
45
46/*
47 * If the previous step in the resolution is reduced, then this check can be
48 * used to determine lower order terms.
49 */
50static inline bool check_variables(const std::vector<bool> &variables,
51        const poly m)
52{
53    const ring R = currRing;
54    // variables[R->N] is true iff index == 1, that is, for the first step in
55    // the resolution
56    if (UNLIKELY(variables[R->N]))
57    {
58        return true;
59    }
60    for (int j = R->N; j > 0; j--)
61    {
62        if (UNLIKELY(variables[j-1] && p_GetExp(m, j, R) > 0))
63        {
64            return true;
65        }
66    }
67    return false;
68}
69
70/*
71 * For each step in the resolution, the following data is saved for each of the
72 * induced leading terms: the leading term itself, its short exponent vector,
73 * and its position in the ideal/module.
74 */
75typedef struct {
76    poly lt;
77    unsigned long sev;
78    unsigned long comp;
79} lt_struct;
80
81static void initialize_hash(lt_struct **C, const ideal L)
82{
83    const ring R = currRing;
84    const unsigned long n_elems = L->ncols;
85    unsigned int *count
86        = (unsigned int *)omAlloc0((L->rank+1)*sizeof(unsigned int));
87    unsigned long k = 0;
88    while (k < n_elems)
89    {
90        count[__p_GetComp(L->m[k], R)]++;
91        k++;
92    }
93    for (int i = 0; i <= L->rank; i++)
94    {
95        // do ++count[i] and use C[i][0].comp to save count[i]
96        C[i] = (lt_struct *)omalloc0((++count[i])*sizeof(lt_struct));
97        C[i][0].comp = count[i];
98    }
99    k = n_elems;
100    // the order of the elements in each C[i] matters if check_variables() is
101    // to be used
102    while (k > 0)
103    {
104        const poly a = L->m[k-1];
105        const unsigned long comp = __p_GetComp(a, R);
106        C[comp][--count[comp]]
107            = (lt_struct){a, p_GetShortExpVector(a, R), k};
108        k--;
109    }
110    omFree(count);
111}
112
113/*
114 * compute a new term in the resolution, that is, compute
115 * ( t * multiplier / f ) where f is an induced leading term from the previous
116 * module, or return NULL if no such f dividing t * multiplier exists, that is,
117 * if multiplier is a lower order term
118 */
119static poly find_reducer(const poly multiplier, const poly t,
120        const lt_struct *const *const hash_previous_module)
121{
122    const ring r = currRing;
123    const lt_struct *v = hash_previous_module[__p_GetComp(t, r)];
124    unsigned long count = v[0].comp;
125    if (UNLIKELY(count == 1))
126    {
127        return NULL;
128    }
129    const poly q = p_New(r);
130    pNext(q) = NULL;
131    p_MemSum_LengthGeneral(q->exp, multiplier->exp, t->exp, r->ExpL_Size);
132    const unsigned long q_not_sev = ~p_GetShortExpVector(q, r);
133    for(unsigned long i = 1; i < count; i++)
134    {
135        if (LIKELY(v[i].sev & q_not_sev)
136                || UNLIKELY(!(_p_LmDivisibleByNoComp(v[i].lt, q, r))))
137                {
138            continue;
139        }
140        p_MemAdd_NegWeightAdjust(q, r);
141        p_ExpVectorDiff(q, q, v[i].lt, r);
142        p_SetComp(q, v[i].comp, r);
143        p_Setm(q, r);
144        number n = n_Div(p_GetCoeff(multiplier, r), p_GetCoeff(v[i].lt, r), r);
145        n_InpMult(n, p_GetCoeff(t, r), r);
146        p_SetCoeff0(q, n_InpNeg(n, r), r);
147        return q;
148    }
149    p_LmFree(q, r);
150    return NULL;
151}
152
153static poly traverse_tail(const poly multiplier, const int comp,
154        const ideal previous_module, const std::vector<bool> &variables,
155        const lt_struct *const *const hash_previous_module);
156
157static poly compute_image(const poly multiplier, const int comp,
158        const ideal previous_module, const std::vector<bool> &variables,
159        const lt_struct *const *const hash_previous_module,
160        const bool use_cache);
161
162/*
163 * recursively call traverse_tail() for each new term found by find_reducer()
164 */
165static poly reduce_term(const poly multiplier, const poly term,
166        const ideal previous_module, const std::vector<bool> &variables,
167        const lt_struct *const *const hash_previous_module,
168        const bool use_cache)
169{
170    poly s = find_reducer(multiplier, term, hash_previous_module);
171    if (s == NULL)
172    {
173        return NULL;
174    }
175    const ring r = currRing;
176    const int c = __p_GetComp(s, r) - 1;
177    poly t;
178    if (use_cache)
179    {
180        t = traverse_tail(s, c, previous_module, variables,
181                hash_previous_module);
182    } else {
183        t = compute_image(s, c, previous_module, variables,
184                hash_previous_module, false);
185    }
186    return p_Add_q(s, t, r);
187}
188
189/*
190 * iterating over tail, call reduce_term(multiplier, p, ...) for each term p in
191 * tail and sum up the results
192 */
193static poly compute_image(const poly multiplier, const int comp,
194        const ideal previous_module, const std::vector<bool> &variables,
195        const lt_struct *const *const hash_previous_module,
196        const bool use_cache)
197{
198    const poly tail = previous_module->m[comp]->next;
199    if (UNLIKELY(tail == NULL) || !check_variables(variables, multiplier))
200    {
201        return NULL;
202    }
203    sBucket_pt sum = sBucketCreate(currRing);
204    for (poly p = tail; p != NULL; p = pNext(p))
205    {
206        const poly rt = reduce_term(multiplier, p, previous_module, variables,
207                hash_previous_module, use_cache);
208        sBucket_Add_p(sum, rt, pLength(rt));
209    }
210    poly s;
211    int l;
212    sBucketClearAdd(sum, &s, &l);
213    sBucketDestroy(&sum);
214    return s;
215}
216
217struct cache_compare
218{
219    inline bool operator() (const poly& l, const poly& r) const
220    {
221        return (p_LmCmp(l, r, currRing) == -1);
222        /* For expensive orderings, consider:
223         * return (memcmp(l->exp, r->exp,
224         *         (currRing->CmpL_Size)*sizeof(unsigned long)) < 0);
225         */
226    }
227};
228
229typedef std::map<poly, poly, cache_compare> cache_term;
230
231static cache_term *Cache;
232
233static void initialize_cache(const int size)
234{
235    Cache = new cache_term[size];
236}
237
238static void delete_cache(const int size)
239{
240    const ring r = currRing;
241    for (int i = 0; i < size; i++)
242    {
243        cache_term *T = &(Cache[i]);
244        for (cache_term::iterator itr = T->begin(); itr != T->end(); ++itr)
245        {
246            p_Delete(&(itr->second), r);
247            p_Delete(const_cast<poly*>(&(itr->first)), r);
248        }
249        T->clear();
250    }
251    delete[](Cache);
252}
253
254static void insert_into_cache_term(cache_term *T, const poly multiplier,
255        const poly p)
256{
257    const ring r = currRing;
258    T->insert(cache_term::value_type(p_Head(multiplier, r), p_Copy(p, r)));
259}
260
261static poly get_from_cache_term(const cache_term::const_iterator itr,
262        const poly multiplier)
263{
264    if (LIKELY(itr->second == NULL))
265    {
266        return NULL;
267    }
268    const ring r = currRing;
269    poly p = p_Copy(itr->second, r);
270    if (LIKELY(!n_Equal(pGetCoeff(multiplier), pGetCoeff(itr->first), r)))
271    {
272        number n = n_Div(pGetCoeff(multiplier), pGetCoeff(itr->first), r);
273        p = p_Mult_nn(p, n, r);
274        n_Delete(&n, r);
275    }
276    return p;
277}
278
279static poly traverse_tail(const poly multiplier, const int comp,
280        const ideal previous_module, const std::vector<bool> &variables,
281        const lt_struct *const *const hash_previous_module)
282{
283    cache_term *T = &(Cache[comp]);
284    cache_term::const_iterator itr = T->find(multiplier);
285    if (LIKELY(itr != T->end()))
286    {
287        return get_from_cache_term(itr, multiplier);
288    }
289    poly p = compute_image(multiplier, comp, previous_module, variables,
290            hash_previous_module, true);
291    insert_into_cache_term(T, multiplier, p);
292    return p;
293}
294
295/*
296 * lift the extended induced leading term a to a syzygy
297 */
298static poly lift_ext_LT(const poly a, const ideal previous_module,
299        const std::vector<bool> &variables,
300        const lt_struct *const *const hash_previous_module,
301        const bool use_cache)
302{
303    const ring R = currRing;
304    // the leading term does not need to be cached
305    poly t1 = compute_image(a, __p_GetComp(a, R)-1, previous_module, variables,
306            hash_previous_module, use_cache);
307    poly t2;
308    if (use_cache)
309    {
310        t2 = traverse_tail(a->next, __p_GetComp(a->next, R)-1,
311                previous_module, variables, hash_previous_module);
312    } else {
313        t2 = compute_image(a->next, __p_GetComp(a->next, R)-1,
314                previous_module, variables, hash_previous_module, false);
315    }
316    t1 = p_Add_q(t1, t2, R);
317    return t1;
318}
319
320/*****************************************************************************/
321
322/*
323 * copied from id_DelDiv(), but without testing and without HAVE_RINGS;
324 * delete id[j], if LT(j) == coeff*mon*LT(i) and vice versa, that is,
325 * delete id[i], if LT(i) == coeff*mon*LT(j)
326 */
327static void id_DelDiv_no_test(ideal id)
328{
329    const ring r = currRing;
330    int i, j;
331    int k = id->ncols-1;
332    for (i = k; i >= 0; i--)
333    {
334        for (j = k; j > i; j--)
335        {
336            if (id->m[j] != NULL)
337            {
338                if (p_DivisibleBy(id->m[i], id->m[j], r))
339                {
340                    p_Delete(&id->m[j], r);
341                }
342                else if (p_DivisibleBy(id->m[j], id->m[i], r))
343                {
344                    p_Delete(&id->m[i], r);
345                    break;
346                }
347            }
348        }
349    }
350}
351
352typedef poly syzHeadFunction(ideal, int, int);
353
354/*
355 * compute the induced leading term corresponding to the index pair (i, j)
356 */
357static poly syzHeadFrame(const ideal G, const int i, const int j)
358{
359    const ring r = currRing;
360    const poly f_i = G->m[i];
361    const poly f_j = G->m[j];
362    poly head = p_Init(r);
363    pSetCoeff0(head, n_Init(1, r->cf));
364    long exp_i, exp_j, lcm;
365    for (int k = (int)r->N; k > 0; k--)
366    {
367        exp_i = p_GetExp(f_i, k, r);
368        exp_j = p_GetExp(f_j, k, r);
369        lcm = si_max(exp_i, exp_j);
370        p_SetExp(head, k, lcm-exp_i, r);
371    }
372    p_SetComp(head, i+1, r);
373    p_Setm(head, r);
374    return head;
375}
376
377/*
378 * compute the _extended_ induced leading term corresponding to the index pair
379 * (i, j), that is, the first two terms w.r.t. the induced order
380 */
381static poly syzHeadExtFrame(const ideal G, const int i, const int j)
382{
383    const ring r = currRing;
384    const poly f_i = G->m[i];
385    const poly f_j = G->m[j];
386    poly head = p_Init(r);
387    pSetCoeff0(head, n_Init(1, r->cf));
388    poly head_ext = p_Init(r);
389    pSetCoeff0(head_ext, n_InpNeg(n_Div(pGetCoeff(f_i), pGetCoeff(f_j), r->cf),
390                r->cf));
391    long exp_i, exp_j, lcm;
392    for (int k = (int)r->N; k > 0; k--)
393    {
394        exp_i = p_GetExp(f_i, k, r);
395        exp_j = p_GetExp(f_j, k, r);
396        lcm = si_max(exp_i, exp_j);
397        p_SetExp(head, k, lcm-exp_i, r);
398        p_SetExp(head_ext, k, lcm-exp_j, r);
399    }
400    p_SetComp(head, i+1, r);
401    p_Setm(head, r);
402    p_SetComp(head_ext, j+1, r);
403    p_Setm(head_ext, r);
404    head->next = head_ext;
405    return head;
406}
407
408typedef ideal syzM_i_Function(ideal, int, syzHeadFunction);
409
410/*
411 * compute the monomial ideal M_i, see reference;
412 * in the first step, we cannot assume that all leading terms which lie in the
413 * component are adjacent to each other
414 */
415static ideal syzM_i_unsorted(const ideal G, const int i,
416        syzHeadFunction *syzHead)
417{
418    const ring r = currRing;
419    ideal M_i = NULL;
420    unsigned long comp = __p_GetComp(G->m[i], r);
421    int ncols = 0;
422    for (int j = i-1; j >= 0; j--)
423    {
424        if (__p_GetComp(G->m[j], r) == comp) ncols++;
425    }
426    if (ncols > 0)
427    {
428        M_i = idInit(ncols, G->ncols);
429        int k = ncols-1;
430        for (int j = i-1; j >= 0; j--)
431        {
432            if (__p_GetComp(G->m[j], r) == comp)
433            {
434                M_i->m[k] = syzHead(G, i, j);
435                k--;
436            }
437        }
438        id_DelDiv_no_test(M_i);
439        idSkipZeroes(M_i);
440    }
441    return M_i;
442}
443
444/*
445 * compute the monomial ideal M_i, see reference;
446 * from step two on, we can assume that all leading terms which lie in the same
447 * component are adjacent to each other
448 */
449static ideal syzM_i_sorted(const ideal G, const int i,
450        syzHeadFunction *syzHead)
451{
452    const ring r = currRing;
453    ideal M_i = NULL;
454    unsigned long comp = __p_GetComp(G->m[i], r);
455    int index = i-1;
456    while (__p_GetComp(G->m[index], r) == comp) index--;
457    index++;
458    int ncols = i-index;
459    if (ncols > 0)
460    {
461        M_i = idInit(ncols, G->ncols);
462        for (int j = ncols-1; j >= 0; j--)
463        {
464            M_i->m[j] = syzHead(G, i, j+index);
465        }
466        id_DelDiv_no_test(M_i);
467        idSkipZeroes(M_i);
468    }
469    return M_i;
470}
471
472/*
473 * concatenate the ideals in M[]
474 */
475static ideal idConcat(const ideal *M, const int size, const int rank)
476{
477    int ncols = 0;
478    for (int i = size-1; i >= 0; i--)
479    {
480        if (M[i] != NULL)
481        {
482            ncols += M[i]->ncols;
483        }
484    }
485    if (ncols == 0) return idInit(1, rank);
486    ideal result = idInit(ncols, rank);
487    int k = ncols-1;
488    for (int i = size-1; i >= 0; i--)
489    {
490        if (M[i] != NULL)
491        {
492            for (int j = M[i]->ncols-1; j >= 0; j--)
493            {
494                result->m[k] = M[i]->m[j];
495                k--;
496            }
497        }
498    }
499    return result;
500}
501
502static int compare_comp(const poly p_a, const poly p_b)
503{
504    const ring r = currRing;
505    long comp_a = __p_GetComp(p_a, r);
506    long comp_b = __p_GetComp(p_b, r);
507    return (comp_a > comp_b) - (comp_a < comp_b);
508}
509
510static int compare_deg(const poly p_a, const poly p_b)
511{
512    const ring r = currRing;
513    long deg_a = p_Deg(p_a, r);
514    long deg_b = p_Deg(p_b, r);
515    return (deg_a > deg_b) - (deg_a < deg_b);
516}
517
518static int compare_lex(const poly p_a, const poly p_b)
519{
520    int cmp;
521    const ring r = currRing;
522    int exp_a[r->N+1];
523    int exp_b[r->N+1];
524    p_GetExpV(p_a, exp_a, r);
525    p_GetExpV(p_b, exp_b, r);
526    for (int i = r->N; i > 0; i--)
527    {
528        cmp = (exp_a[i] > exp_b[i]) - (exp_a[i] < exp_b[i]);
529        if (cmp != 0)
530        {
531            return cmp;
532        }
533    }
534    return 0;
535}
536
537static int compare_Mi(const void* a, const void *b)
538{
539    poly p_a = *((poly *)a);
540    poly p_b = *((poly *)b);
541    int cmp;
542    if ((cmp = compare_comp(p_a, p_b))
543            || (cmp = compare_deg(p_a, p_b))
544            || (cmp = compare_lex(p_a, p_b)))
545            {
546        return cmp;
547    }
548    return 0;
549}
550
551/*
552 * compute the frame, that is, the induced leading terms for the next step in
553 * the resolution
554 */
555static ideal computeFrame(const ideal G, syzM_i_Function syzM_i,
556        syzHeadFunction *syzHead)
557{
558    ideal *M = (ideal *)omalloc((G->ncols-1)*sizeof(ideal));
559    for (int i = G->ncols-2; i >= 0; i--)
560    {
561        M[i] = syzM_i(G, i+1, syzHead);
562    }
563    ideal frame = idConcat(M, G->ncols-1, G->ncols);
564    for (int i = G->ncols-2; i >= 0; i--)
565    {
566        if (M[i] != NULL)
567        {
568            omFreeSize(M[i]->m, M[i]->ncols*sizeof(poly));
569            omFreeBin(M[i], sip_sideal_bin);
570        }
571    }
572    omfree(M);
573    qsort(frame->m, frame->ncols, sizeof(poly), compare_Mi);
574    return frame;
575}
576
577/*
578 * lift each (extended) induced leading term to a syzygy
579 */
580static void computeLiftings(const resolvente res, const int index,
581        const std::vector<bool> &variables, const bool use_cache)
582{
583    if (use_cache)
584    {
585        initialize_cache(res[index-1]->ncols);
586    }
587    lt_struct **hash_previous_module
588        = (lt_struct **)omAlloc((res[index-1]->rank+1)*sizeof(lt_struct *));
589    initialize_hash(hash_previous_module, res[index-1]);
590    for (int j = res[index]->ncols-1; j >= 0; j--)
591    {
592        res[index]->m[j]->next->next = lift_ext_LT(res[index]->m[j],
593                res[index-1], variables, hash_previous_module, use_cache);
594    }
595    for (int i = 0; i <= res[index-1]->rank; i++)
596    {
597        omfree(hash_previous_module[i]);
598    }
599    omFree(hash_previous_module);
600    if (use_cache)
601    {
602        delete_cache(res[index-1]->ncols);
603    }
604}
605
606/*
607 * check if the monomial m contains any of the variables set to false
608 */
609static inline bool contains_unused_variable(const poly m,
610    const std::vector<bool> &variables)
611{
612    const ring R = currRing;
613    for (int j = R->N; j > 0; j--)
614    {
615        if (!variables[j-1] && p_GetExp(m, j, R) > 0)
616        {
617            return true;
618        }
619    }
620    return false;
621}
622
623/*
624 * delete any term in res[index] which contains any of the variables set to
625 * false
626 */
627static void delete_variables(resolvente res, const int index,
628    const std::vector<bool> &variables)
629{
630    for (int i = 0; i < res[index]->ncols; i++)
631    {
632        poly p_iter = res[index]->m[i]->next;
633        if (p_iter != NULL)
634        {
635            while (p_iter->next != NULL)
636            {
637                if (contains_unused_variable(p_iter->next, variables))
638                {
639                    pLmDelete(&p_iter->next);
640                } else {
641                    pIter(p_iter);
642                }
643            }
644        }
645    }
646}
647
648static void delete_tails(resolvente res, const int index)
649{
650    const ring r = currRing;
651    for (int i = 0; i < res[index]->ncols; i++)
652    {
653        if (res[index]->m[i] != NULL)
654        {
655            p_Delete(&(res[index]->m[i]->next), r);
656        }
657    }
658}
659
660/*
661 * for each step in the resolution, compute the corresponding module until
662 * either index == max_index is reached or res[index] is the zero module
663 */
664static int computeResolution_iteration(resolvente res, const int max_index,
665        syzHeadFunction *syzHead, const bool do_lifting,
666        const bool single_module, const bool use_cache,
667        const bool use_tensor_trick, std::vector<bool> &variables)
668{
669    int index = 1;
670    while (!idIs0(res[index]))
671    {
672        if (do_lifting)
673        {
674            computeLiftings(res, index, variables, use_cache);
675            if (single_module)
676            {
677                delete_tails(res, index-1);
678            }
679            // we don't know if the input is a reduced SB:
680            if (index == 1)
681            {
682                variables[currRing->N] = false;
683            }
684            update_variables(variables, res[index]);
685            if (use_tensor_trick)
686            {
687                delete_variables(res, index, variables);
688            }
689        }
690        if (index >= max_index) { break; }
691        index++;
692        res[index] = computeFrame(res[index-1], syzM_i_sorted, syzHead);
693    }
694    return index;
695}
696
697/*
698 * compute the frame of the first syzygy module and set variables, then call
699 * computeResolution_iteration() for the remaining steps
700 */
701static int computeResolution(resolvente res, const int max_index,
702        syzHeadFunction *syzHead, const bool do_lifting,
703        const bool single_module, const bool use_cache,
704        const bool use_tensor_trick)
705{
706    if (idIs0(res[0]))
707    {
708        return 1;
709    }
710    std::vector<bool> variables;
711    variables.resize(currRing->N+1, true);
712    if (do_lifting)
713    {
714        update_variables(variables, res[0]);
715        if (use_tensor_trick)
716        {
717            delete_variables(res, 0, variables);
718        }
719    }
720    int index = 0;
721    if (max_index > 0)
722    {
723        res[1] = computeFrame(res[0], syzM_i_unsorted, syzHead);
724        index = computeResolution_iteration(res, max_index, syzHead,
725                do_lifting, single_module, use_cache, use_tensor_trick,
726                variables);
727    }
728    variables.clear();
729    return index+1;
730}
731
732static void set_options(syzHeadFunction **syzHead_ptr, bool *do_lifting_ptr,
733        bool *single_module_ptr, const char *method)
734{
735    if (strcmp(method, "complete") == 0)
736    {   // default
737        *syzHead_ptr = syzHeadExtFrame;
738        *do_lifting_ptr = true;
739        *single_module_ptr = false;
740    }
741    else if (strcmp(method, "frame") == 0)
742    {
743        *syzHead_ptr = syzHeadFrame;
744        *do_lifting_ptr = false;
745        *single_module_ptr = false;
746    }
747    else if (strcmp(method, "extended frame") == 0)
748    {
749        *syzHead_ptr = syzHeadExtFrame;
750        *do_lifting_ptr = false;
751        *single_module_ptr = false;
752    }
753    else if (strcmp(method, "single module") == 0)
754    {
755        *syzHead_ptr = syzHeadExtFrame;
756        *do_lifting_ptr = true;
757        *single_module_ptr = true;
758    }
759    else {   // "linear strand" (not yet implemented)
760        *syzHead_ptr = syzHeadExtFrame;
761        *do_lifting_ptr = true;
762        *single_module_ptr = false;
763    }
764}
765
766/*
767 * insert the first term of r at the right place
768 */
769#define insert_first_term(r, p, q, R)                             \
770do                                                                \
771{                                                                 \
772    p = r;                                                        \
773    q = p->next;                                                  \
774    if (q != NULL && p_LmCmp(p, q, R) != 1) {                     \
775        while (q->next != NULL && p_LmCmp(p, q->next, R) == -1) { \
776            pIter(q);                                             \
777        }                                                         \
778        r = p->next;                                              \
779        p->next = q->next;                                        \
780        q->next = p;                                              \
781    }                                                             \
782}                                                                 \
783while (0)
784
785/*
786 * For each poly in the resolution, insert the first two terms at their right
787 * places. If single_module is true, then only consider the last module.
788 */
789static void insert_ext_induced_LTs(const resolvente res, const int length,
790        const bool single_module)
791{
792    const ring R = currRing;
793    poly p, q;
794    int index = (single_module ? length-1 : 1);
795    while (index < length && !idIs0(res[index]))
796    {
797        for (int j = res[index]->ncols-1; j >= 0; j--)
798        {
799            insert_first_term(res[index]->m[j]->next, p, q, R);
800            insert_first_term(res[index]->m[j], p, q, R);
801        }
802        index++;
803    }
804}
805
806/*
807 * Compute the Schreyer resolution of arg, see reference at the beginning of
808 * this file.
809 *
810 * If use_cache == true (default), the result of compute_image() is cached for
811 * _every_ term in the current step of the resolution. This corresponds to the
812 * subtree attached to the node which represents this term, see reference.
813 *
814 * If use_tensor_trick == true, the current module is modfied after each
815 * lifting step in the resolution: any term which contains a variable which
816 * does not appear among the (induced) leading terms is deleted. Note that the
817 * resulting object is not necessarily a complex anymore. However, constant
818 * entries remain exactly the same. This option does not apply for
819 * method == "frame" and method "extended frame".
820 *
821 * These two options are used in PrymGreen.jl; do not delete!
822 */
823syStrategy syFrank(const ideal arg, const int length, const char *method,
824        const bool use_cache, const bool use_tensor_trick)
825{
826    syStrategy result = (syStrategy)omAlloc0(sizeof(ssyStrategy));
827    resolvente res = (resolvente)omAlloc0((length+1)*sizeof(ideal));
828    if (strcmp(method, "frame") != 0)
829    {
830        res[0] = id_Copy(arg, currRing);
831    }
832    else
833    {
834        res[0] = id_Head(arg, currRing);
835    }
836    syzHeadFunction *syzHead;
837    bool do_lifting;
838    bool single_module;
839    set_options(&syzHead, &do_lifting, &single_module, method);
840    int new_length = computeResolution(res, length-1, syzHead, do_lifting,
841            single_module, use_cache, use_tensor_trick);
842    if (new_length < length)
843    {
844        res = (resolvente)omReallocSize(res, (length+1)*sizeof(ideal),
845                (new_length+1)*sizeof(ideal));
846    }
847    if (strcmp(method, "frame") != 0)
848    {
849        insert_ext_induced_LTs(res, new_length, single_module);
850    }
851    result->fullres = res;
852    result->length = new_length;
853    result->list_length = new_length;
854    return result;
855}
856
Note: See TracBrowser for help on using the repository browser.