source: git/kernel/GBEngine/syz4.cc @ 4c6e9d4

fieker-DuValspielwiese
Last change on this file since 4c6e9d4 was 4c6e9d4, checked in by Andreas Steenpass <steenpass@…>, 8 years ago
chg: replace TCacheKey_test and TCacheValue_test by poly
  • Property mode set to 100644
File size: 20.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: resolutions
6*/
7
8#include <kernel/GBEngine/syz.h>
9#include <omalloc/omalloc.h>
10#include <coeffs/numbers.h>
11#include <kernel/polys.h>
12#include <kernel/ideals.h>
13
14#include <vector>
15#include <map>
16
17typedef struct {
18    poly lt;
19    unsigned long sev;
20    unsigned int label;
21} CLeadingTerm_struct;
22
23typedef std::vector<const CLeadingTerm_struct*> TReducers_test;
24typedef std::map<long, TReducers_test> CReducersHash_test;
25
26static poly TraverseTail_test(poly multiplier, const int tail,
27   const ideal m_idTails_test, const std::vector<bool> &m_lcm,
28   const CReducersHash_test *m_div, const CReducersHash_test *m_checker);
29static poly ComputeImage_test(poly multiplier, const int tail,
30    const ideal m_idTails_test, const std::vector<bool> &m_lcm,
31    const CReducersHash_test *m_div, const CReducersHash_test *m_checker);
32static poly TraverseTail_test(poly multiplier, poly tail,
33    const ideal m_idTails_test, const std::vector<bool> &m_lcm,
34    const CReducersHash_test *m_div, const CReducersHash_test *m_checker);
35static inline poly ReduceTerm_test(poly multiplier, poly term4reduction,
36    poly syztermCheck, const ideal m_idTails_test,
37    const std::vector<bool> &m_lcm, const CReducersHash_test *m_div,
38    const CReducersHash_test *m_checker);
39static poly leadmonom_test(const poly p, const ring r, const bool bSetZeroComp = true);
40
41static std::vector<bool> CLCM_test_redefine(const ideal L)
42{
43  std::vector<bool> clcm;
44  const ring R = currRing;
45  if( L != NULL )
46  {
47    const int l = IDELEMS(L);
48    clcm.resize(currRing->N, false);
49    for( int k = l - 1; k >= 0; k-- )
50    {
51      const poly a = L->m[k];
52      for (unsigned int j = currRing->N; j > 0; j--)
53        if ( !clcm[j-1] )
54          clcm[j-1] = (p_GetExp(a, j, R) > 0);
55    }
56  }
57  return clcm;
58}
59
60bool CLCM_test_Check(const std::vector<bool> &clcm, const poly m)
61{
62  if(m != NULL)
63  {
64    const ring R = currRing;
65    for (unsigned int j = currRing->N; j > 0; j--)
66      if ( clcm[j-1] )
67        if(p_GetExp(m, j, R) > 0)
68          return true;
69    return false;
70  } else return true;
71}
72
73static poly TraverseNF_test(const poly a, ideal m_idLeads_test,
74    const ideal m_idTails_test, const std::vector<bool> &m_lcm,
75    const CReducersHash_test *m_div, const CReducersHash_test *m_checker)
76{
77  const ideal& L = m_idLeads_test;
78  const ring R = currRing;
79  const int r = p_GetComp(a, R) - 1;
80  poly aa = leadmonom_test(a, R);
81  poly t = TraverseTail_test(aa, r, m_idTails_test, m_lcm, m_div, m_checker);
82  t = p_Add_q(t,
83      ReduceTerm_test(aa, L->m[r], a, m_idTails_test, m_lcm, m_div, m_checker),
84      R);
85  p_Delete(&aa, R);
86  return t;
87}
88
89#define CACHE 1
90
91#if CACHE
92bool my_p_LmCmp_test (poly a, poly b, const ring r)
93{
94  return p_LmCmp(a, b, r) == -1;
95}
96
97struct CCacheCompare_test
98{
99  const ring& m_ring_test;
100  CCacheCompare_test(): m_ring_test(currRing) {}
101  CCacheCompare_test(const ring& r): m_ring_test(r) {}
102  CCacheCompare_test(const CCacheCompare_test& lhs):
103    m_ring_test(lhs.m_ring_test) {}
104  CCacheCompare_test& operator=(const CCacheCompare_test& lhs)
105  {
106    return (const_cast<CCacheCompare_test&>(lhs));
107  }
108  inline bool operator() (const poly& l, const poly& r)
109    const
110  {
111    return my_p_LmCmp_test(l, r, m_ring_test);
112  }
113};
114
115typedef std::map<poly, poly, CCacheCompare_test>
116  TP2PCache_test;
117typedef std::map<int, TP2PCache_test> TCache_test;
118
119static TCache_test m_cache_test;
120
121static FORCE_INLINE poly myp_Head_test(const poly p, const bool bIgnoreCoeff,
122  const ring r)
123{
124  p_LmCheckPolyRing1(p, r);
125  poly np;
126  omTypeAllocBin(poly, np, r->PolyBin);
127  p_SetRingOfLm(np, r);
128  memcpy(np->exp, p->exp, r->ExpL_Size*sizeof(long));
129  pNext(np) = NULL;
130  pSetCoeff0(np, (bIgnoreCoeff)? NULL : n_Copy(pGetCoeff(p), r->cf));
131  p_LmCheckPolyRing1(np, r);
132  return np;
133}
134#endif   // CACHE
135
136static poly TraverseTail_test(poly multiplier, const int tail,
137    const ideal m_idTails_test, const std::vector<bool> &m_lcm,
138    const CReducersHash_test *m_div, const CReducersHash_test *m_checker)
139{
140  const ring& r = currRing;
141#if CACHE
142  TCache_test::iterator top_itr = m_cache_test.find(tail);
143  if ( top_itr != m_cache_test.end() )
144  {
145     TP2PCache_test& T = top_itr->second;
146     TP2PCache_test::iterator itr = T.find(multiplier);
147     if( itr != T.end() )
148     {
149       if( itr->second == NULL )
150         return (NULL);
151       poly p = p_Copy(itr->second, r);
152       if( !n_Equal( pGetCoeff(multiplier), pGetCoeff(itr->first), r) )
153       {
154         number n = n_Div( pGetCoeff(multiplier), pGetCoeff(itr->first), r);
155         p = p_Mult_nn(p, n, r);
156         n_Delete(&n, r);
157       }
158       return p;
159     }
160     const poly p = ComputeImage_test(multiplier, tail, m_idTails_test, m_lcm,
161         m_div, m_checker);
162     itr = T.find(multiplier);
163     if( itr == T.end() )
164     {
165       T.insert(TP2PCache_test::value_type(myp_Head_test(multiplier, (p==NULL),
166         r), p) );
167       return p_Copy(p, r);
168     }
169     return p;
170  }
171#endif   // CACHE
172  const poly p = ComputeImage_test(multiplier, tail, m_idTails_test, m_lcm,
173      m_div, m_checker);
174#if CACHE
175  top_itr = m_cache_test.find(tail);
176  if ( top_itr != m_cache_test.end() )
177  {
178    TP2PCache_test& T = top_itr->second;
179    TP2PCache_test::iterator itr = T.find(multiplier);
180    if( itr == T.end() )
181    {
182      T.insert(TP2PCache_test::value_type(myp_Head_test(multiplier, (p==NULL),
183          r), p));
184      return p_Copy(p, r);
185    }
186    return p;
187  }
188  CCacheCompare_test o(r);
189  TP2PCache_test T(o);
190  T.insert(TP2PCache_test::value_type(myp_Head_test(multiplier, (p==NULL), r),
191    p));
192  m_cache_test.insert( TCache_test::value_type(tail, T) );
193  return p_Copy(p, r);
194#else
195  return p;
196#endif   // CACHE
197}
198
199#define BUCKETS 1
200/* 0: use original code (not implemented)
201 * 1: use Singular's SBuckets
202 * 2: use Singular's polynomial arithmetic
203 */
204
205static poly ComputeImage_test(poly multiplier, const int t,
206    const ideal m_idTails_test, const std::vector<bool> &m_lcm,
207    const CReducersHash_test *m_div, const CReducersHash_test *m_checker)
208{
209  const poly tail = m_idTails_test->m[t];
210  if(tail != NULL)
211  {
212    const ring r = currRing;
213    if( !CLCM_test_Check(m_lcm, multiplier) )
214    {
215      return NULL;
216    }
217#if BUCKETS == 0
218    SBucketWrapper sum(r, m_sum_bucket_factory);
219#elif BUCKETS == 1
220    sBucket_pt sum = sBucketCreate(currRing);
221#else   // BUCKETS == 2
222    poly s = NULL;
223#endif   // BUCKETS
224    for(poly p = tail; p != NULL; p = pNext(p))   // iterate over the tail
225    {
226      const poly rt = ReduceTerm_test(multiplier, p, NULL, m_idTails_test,
227          m_lcm, m_div, m_checker);
228#if BUCKETS == 0
229      sum.Add(rt);
230#elif BUCKETS == 1
231      sBucket_Add_p(sum, rt, pLength(rt));
232#else   // BUCKETS == 2
233      s = p_Add_q(s, rt, r);
234#endif   // BUCKETS
235    }
236#if BUCKETS == 0
237    const poly s = sum.ClearAdd();
238#elif BUCKETS == 1
239    poly s;
240    int l;
241    sBucketClearAdd(sum, &s, &l);
242    sBucketDestroy(&sum);
243#else   // BUCKETS == 2
244    // nothing to do
245#endif   // BUCKETS
246    return s;
247  }
248  return NULL;
249}
250
251static poly leadmonom_test(const poly p, const ring r, const bool bSetZeroComp)
252{
253  if( p == NULL )
254    return NULL;
255  poly m = p_LmInit(p, r);
256  p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);
257  if( bSetZeroComp )
258    p_SetComp(m, 0, r);
259  p_Setm(m, r);
260  return m;
261}
262
263// _p_LmDivisibleByNoComp for a | b*c
264static inline BOOLEAN _p_LmDivisibleByNoComp(const poly a, const poly b, const poly c, const ring r)
265{
266  int i=r->VarL_Size - 1;
267  unsigned long divmask = r->divmask;
268  unsigned long la, lb;
269
270  if (r->VarL_LowIndex >= 0)
271  {
272    i += r->VarL_LowIndex;
273    do
274    {
275      la = a->exp[i];
276      lb = b->exp[i] + c->exp[i];
277      if ((la > lb) ||
278          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
279      {
280        return FALSE;
281      }
282      i--;
283    }
284    while (i>=r->VarL_LowIndex);
285  }
286  else
287  {
288    do
289    {
290      la = a->exp[r->VarL_Offset[i]];
291      lb = b->exp[r->VarL_Offset[i]] + c->exp[r->VarL_Offset[i]];
292      if ((la > lb) ||
293          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
294      {
295        return FALSE;
296      }
297      i--;
298    }
299    while (i>=0);
300  }
301  return TRUE;
302}
303
304static void deleteCRH(CReducersHash_test *C)
305{
306    for (CReducersHash_test::iterator it = C->begin(); it != C->end(); it++) {
307        TReducers_test& v = it->second;
308        for (TReducers_test::const_iterator vit = v.begin(); vit != v.end();
309            vit++) {
310            omfree(const_cast<CLeadingTerm_struct*>(*vit));
311        }
312        v.erase(v.begin(), v.end());
313    }
314    C->erase(C->begin(), C->end());
315}
316
317bool IsDivisible(const CReducersHash_test *C, const poly product)
318{
319    CReducersHash_test::const_iterator m_itr
320        = C->find(p_GetComp(product, currRing));
321    if (m_itr == C->end()) {
322        return false;
323    }
324    TReducers_test::const_iterator m_current = (m_itr->second).begin();
325    TReducers_test::const_iterator m_finish  = (m_itr->second).end();
326    const unsigned long m_not_sev = ~p_GetShortExpVector(product, currRing);
327    for ( ; m_current != m_finish; ++m_current) {
328        if (p_LmShortDivisibleByNoComp((*m_current)->lt, (*m_current)->sev,
329            product, m_not_sev, currRing)) {
330            return true;
331        }
332    }
333    return false;
334}
335
336
337poly FindReducer(const poly multiplier, const poly t, const poly syzterm,
338    const CReducersHash_test *syz_checker, const CReducersHash_test *m_div)
339{
340  const ring r = currRing;
341  CReducersHash_test::const_iterator m_itr
342      = m_div->find(p_GetComp(t, currRing));
343  if (m_itr == m_div->end()) {
344    return NULL;
345  }
346  TReducers_test::const_iterator m_current = (m_itr->second).begin();
347  TReducers_test::const_iterator m_finish  = (m_itr->second).end();
348  if (m_current == m_finish) {
349    return NULL;
350  }
351  const poly q = p_New(r);
352  pNext(q) = NULL;
353  const unsigned long m_not_sev = ~p_GetShortExpVector(multiplier, t, r);
354  for( ; m_current != m_finish; ++m_current) {
355    if ( ((*m_current)->sev & m_not_sev)
356        || !(_p_LmDivisibleByNoComp((*m_current)->lt, multiplier, t, r))) {
357      continue;
358    }
359    const poly p = (*m_current)->lt;
360    const int k  = (*m_current)->label;
361    p_ExpVectorSum(q, multiplier, t, r); // q == product == multiplier * t
362    p_ExpVectorDiff(q, q, p, r); // (LM(product) / LM(L[k]))
363    p_SetComp(q, k + 1, r);
364    p_Setm(q, r);
365    // cannot allow something like: a*gen(i) - a*gen(i)
366    if (syzterm != NULL && (k == p_GetComp(syzterm, r) - 1)
367        && p_ExpVectorEqual(syzterm, q, r)) {
368      continue;
369    }
370    if (IsDivisible(syz_checker, q)) {
371      continue;
372    }
373    number n = n_Mult(p_GetCoeff(multiplier, r), p_GetCoeff(t, r), r);
374    p_SetCoeff0(q, n_InpNeg(n, r), r);
375    return q;
376  }
377  p_LmFree(q, r);
378  return NULL;
379}
380
381static inline poly ReduceTerm_test(poly multiplier, poly term4reduction,
382    poly syztermCheck, const ideal m_idTails_test,
383    const std::vector<bool> &m_lcm, const CReducersHash_test *m_div,
384    const CReducersHash_test *m_checker)
385{
386  const ring r = currRing;
387  poly s = NULL;
388  if( CLCM_test_Check(m_lcm, multiplier) )
389  {
390    s = FindReducer(multiplier, term4reduction, syztermCheck, m_checker,
391        m_div);
392  }
393  if( s == NULL )
394  {
395    return NULL;
396  }
397  poly b = leadmonom_test(s, r);
398  const int c = p_GetComp(s, r) - 1;
399  const poly t
400      = TraverseTail_test(b, c, m_idTails_test, m_lcm, m_div, m_checker);
401  pDelete(&b);
402  if( t != NULL )
403    s = p_Add_q(s, t, r);
404  return s;
405}
406
407static void initialize(CReducersHash_test &C, const ideal L)
408{
409  if( L != NULL )
410  {
411    const ring R = currRing;
412    for( int k = IDELEMS(L) - 1; k >= 0; k-- )
413    {
414      const poly a = L->m[k];
415      if( a != NULL )
416      {
417        CLeadingTerm_struct *CLT
418            = (CLeadingTerm_struct*)omalloc(sizeof(CLeadingTerm_struct));
419        CLT->lt = a;
420        CLT->sev = p_GetShortExpVector(a, R);
421        CLT->label = k;
422        C[p_GetComp(a, R)].push_back( CLT );
423      }
424    }
425  }
426}
427
428/*****************************************************************************/
429
430typedef poly syzHeadFunction(ideal, int, int);
431
432static poly syzHeadFrame(const ideal G, const int i, const int j)
433{
434    const ring r = currRing;
435    const poly f_i = G->m[i];
436    const poly f_j = G->m[j];
437    poly head = p_Init(r);
438    pSetCoeff0(head, n_Init(1, r->cf));
439    long exp_i, exp_j, lcm;
440    for (int k = (int)r->N; k > 0; k--) {
441        exp_i = p_GetExp(f_i, k, r);
442        exp_j = p_GetExp(f_j, k, r);
443        lcm = si_max(exp_i, exp_j);
444        p_SetExp(head, k, lcm-exp_i, r);
445    }
446    p_SetComp(head, i+1, r);
447    p_Setm(head, r);
448    return head;
449}
450
451static poly syzHeadExtFrame(const ideal G, const int i, const int j)
452{
453    const ring r = currRing;
454    const poly f_i = G->m[i];
455    const poly f_j = G->m[j];
456    poly head = p_Init(r);
457    pSetCoeff0(head, n_Init(1, r->cf));
458    poly head_ext = p_Init(r);
459    pSetCoeff0(head_ext, n_InpNeg(n_Div(pGetCoeff(f_j), pGetCoeff(f_i), r->cf),
460        r->cf));
461    long exp_i, exp_j, lcm;
462    for (int k = (int)r->N; k > 0; k--) {
463        exp_i = p_GetExp(f_i, k, r);
464        exp_j = p_GetExp(f_j, k, r);
465        lcm = si_max(exp_i, exp_j);
466        p_SetExp(head, k, lcm-exp_i, r);
467        p_SetExp(head_ext, k, lcm-exp_j, r);
468    }
469    p_SetComp(head, i+1, r);
470    p_Setm(head, r);
471    p_SetComp(head_ext, j+1, r);
472    p_Setm(head_ext, r);
473    head->next = head_ext;
474    return head;
475}
476
477typedef ideal syzM_i_Function(ideal, int, syzHeadFunction);
478
479static ideal syzM_i_unsorted(const ideal G, const int i,
480    const syzHeadFunction *syzHead)
481{
482    ideal M_i = NULL;
483    int comp = pGetComp(G->m[i]);
484    int ncols = 0;
485    for (int j = i-1; j >= 0; j--) {
486        if (pGetComp(G->m[j]) == comp) ncols++;
487    }
488    if (ncols > 0) {
489        M_i = idInit(ncols, G->ncols);
490        int k = ncols-1;
491        for (int j = i-1; j >= 0; j--) {
492            if (pGetComp(G->m[j]) == comp) {
493                M_i->m[k] = syzHead(G, i, j);
494                k--;
495            }
496        }
497        id_DelDiv(M_i, currRing);
498        idSkipZeroes(M_i);
499    }
500    return M_i;
501}
502
503static ideal syzM_i_sorted(const ideal G, const int i,
504    const syzHeadFunction *syzHead)
505{
506    ideal M_i = NULL;
507    int comp = pGetComp(G->m[i]);
508    int index = i-1;
509    while (pGetComp(G->m[index]) == comp) index--;
510    index++;
511    int ncols = i-index;
512    if (ncols > 0) {
513        M_i = idInit(ncols, G->ncols);
514        for (int j = ncols-1; j >= 0; j--) {
515            M_i->m[j] = syzHead(G, i, j+index);
516        }
517        id_DelDiv(M_i, currRing);
518        idSkipZeroes(M_i);
519    }
520    return M_i;
521}
522
523static ideal idConcat(const ideal* M, const int size, const int rank)
524{
525    int ncols = 0;
526    for (int i = size-1; i >= 0; i--) {
527        if (M[i] != NULL) {
528            ncols += M[i]->ncols;
529        }
530    }
531    if (ncols == 0) return idInit(1, rank);
532    ideal result = idInit(ncols, rank);
533    int k = ncols-1;
534    for (int i = size-1; i >= 0; i--) {
535        if (M[i] != NULL) {
536            for (int j = M[i]->ncols-1; j >= 0; j--) {
537                result->m[k] = M[i]->m[j];
538                k--;
539            }
540        }
541    }
542    return result;
543}
544
545#define SORT_MI 1
546
547#if SORT_MI
548static int compare_Mi(const void* a, const void *b)
549{
550    const ring r = currRing;
551    poly p_a = *((poly *)a);
552    poly p_b = *((poly *)b);
553    int cmp;
554    int deg_a = p_Deg(p_a, r);
555    int deg_b = p_Deg(p_b, r);
556    cmp = (deg_a > deg_b) - (deg_a < deg_b);
557    if (cmp != 0) {
558         return cmp;
559    }
560    int comp_a = p_GetComp(p_a, r);
561    int comp_b = p_GetComp(p_b, r);
562    cmp = (comp_a > comp_b) - (comp_a < comp_b);
563    if (cmp != 0) {
564         return cmp;
565    }
566    int exp_a[r->N+1];
567    int exp_b[r->N+1];
568    p_GetExpV(p_a, exp_a, r);
569    p_GetExpV(p_b, exp_b, r);
570    for (int i = r->N; i > 0; i--) {
571        cmp = (exp_a[i] > exp_b[i]) - (exp_a[i] < exp_b[i]);
572        if (cmp != 0) {
573            return cmp;
574        }
575    }
576    return 0;
577}
578#endif   // SORT_MI
579
580static ideal computeFrame(const ideal G, const syzM_i_Function syzM_i,
581    const syzHeadFunction *syzHead)
582{
583    ideal* M = (ideal *)omalloc((G->ncols-1)*sizeof(ideal));
584    for (int i = G->ncols-2; i >= 0; i--) {
585        M[i] = syzM_i(G, i+1, syzHead);
586    }
587    ideal frame = idConcat(M, G->ncols-1, G->ncols);
588    for (int i = G->ncols-2; i >= 0; i--) {
589        if (M[i] != NULL) {
590            omFreeSize(M[i]->m, IDELEMS(M[i])*sizeof(poly));
591            omFreeBin(M[i], sip_sideal_bin);
592        }
593    }
594    omFree(M);
595#if SORT_MI
596    qsort(frame->m, IDELEMS(frame), sizeof(poly), compare_Mi);
597#endif
598    return frame;
599}
600
601static void setGlobalVariables(const resolvente res, const int index,
602    const ideal m_idLeads_test, const ideal m_syzLeads_test)
603{
604#if CACHE
605    for (TCache_test::iterator it = m_cache_test.begin();
606        it != m_cache_test.end(); it++) {
607        TP2PCache_test& T = it->second;
608        for (TP2PCache_test::iterator vit = T.begin(); vit != T.end(); vit++) {
609            p_Delete((&(vit->second)), currRing);
610            p_Delete(const_cast<poly*>(&(vit->first)), currRing);
611        }
612        T.erase(T.begin(), T.end());
613    }
614    m_cache_test.erase(m_cache_test.begin(), m_cache_test.end());
615#endif   // CACHE
616}
617
618#define MYLIFT 0
619
620static void computeLiftings(const resolvente res, const int index)
621{
622#if MYLIFT
623    for (int j = res[index]->ncols-1; j >= 0; j--) {
624        liftTree(res[index]->m[j], res[index-1]);
625    }
626#else
627    ideal m_idLeads_test = idCopy(res[index-1]);
628    ideal m_idTails_test = idInit(IDELEMS(res[index-1]), res[index-1]->rank);
629    for (int i = IDELEMS(res[index-1])-1; i >= 0; i--) {
630        m_idTails_test->m[i] = m_idLeads_test->m[i]->next;
631        m_idLeads_test->m[i]->next = NULL;
632    }
633    std::vector<bool> m_lcm = CLCM_test_redefine(m_idLeads_test);
634    CReducersHash_test m_div;
635    initialize(m_div, m_idLeads_test);
636    ideal m_syzLeads_test = idCopy(res[index]);
637    for (int i = IDELEMS(res[index])-1; i >= 0; i--) {
638        pDelete(&m_syzLeads_test->m[i]->next);
639        m_syzLeads_test->m[i]->next = NULL;
640    }
641    CReducersHash_test m_checker;
642    initialize(m_checker, m_syzLeads_test);
643    setGlobalVariables(res, index, m_idLeads_test, m_syzLeads_test);
644    poly p;
645    for (int j = res[index]->ncols-1; j >= 0; j--) {
646        p = res[index]->m[j];
647        pDelete(&res[index]->m[j]->next);
648        p->next = NULL;
649        res[index]->m[j]->next = TraverseNF_test(p, m_idLeads_test,
650            m_idTails_test, m_lcm, &m_div, &m_checker);
651    }
652    deleteCRH(&m_checker);
653    deleteCRH(&m_div);
654    m_lcm.clear();
655    idDelete(&m_idLeads_test);
656    idDelete(&m_idTails_test);
657    idDelete(&m_syzLeads_test);
658#endif   // MYLIFT
659}
660
661static void sortPolysTails(const resolvente res, const int index)
662{
663    const ring r = currRing;
664    for (int j = res[index]->ncols-1; j >= 0; j--) {
665        if (res[index]->m[j]->next != NULL) {
666            res[index]->m[j]->next->next
667                = p_SortAdd(res[index]->m[j]->next->next, r);
668        }
669    }
670}
671
672static int computeResolution(resolvente &res, const int length,
673    const syzHeadFunction *syzHead)
674{
675    int index = 0;
676    if (!idIs0(res[index]) && index < length) {
677        index++;
678        res[index] = computeFrame(res[index-1], syzM_i_unsorted, syzHead);
679    }
680    while (!idIs0(res[index]) && index < length) {
681#if 1
682        computeLiftings(res, index);
683        sortPolysTails(res, index);
684#endif   // LIFT
685        index++;
686        res[index] = computeFrame(res[index-1], syzM_i_sorted, syzHead);
687    }
688    if (index < length) {
689        res = (resolvente)omReallocSize(res, (length+1)*sizeof(ideal),
690            (index+1)*sizeof(ideal));
691    }
692    return index;
693}
694
695static void sortPolys(const resolvente res, const int length)
696{
697    const ring r = currRing;
698    for (int i = length; i > 0; i--) {
699        for (int j = res[i]->ncols-1; j >= 0; j--) {
700            res[i]->m[j] = p_SortAdd(res[i]->m[j], r, TRUE);
701        }
702    }
703}
704
705syStrategy syFrank(const ideal arg, int &length, const char *method)
706{
707    syStrategy result = (syStrategy)omAlloc0(sizeof(ssyStrategy));
708    resolvente res = (resolvente)omAlloc0((length+1)*sizeof(ideal));
709    res[0] = id_Copy(arg, currRing);
710    syzHeadFunction *syzHead;
711    if (strcmp(method, "frame") == 0 || strcmp(method, "linear strand") == 0) {
712        syzHead = syzHeadFrame;
713    }
714    else {
715        syzHead = syzHeadExtFrame;
716    }
717    length = computeResolution(res, length, syzHead);
718    sortPolys(res, length);
719#if CACHE
720    for (TCache_test::iterator it = m_cache_test.begin();
721        it != m_cache_test.end(); it++) {
722        TP2PCache_test& T = it->second;
723        for (TP2PCache_test::iterator vit = T.begin(); vit != T.end(); vit++) {
724            p_Delete((&(vit->second)), currRing);
725            p_Delete(const_cast<poly*>(&(vit->first)), currRing);
726        }
727        T.erase(T.begin(), T.end());
728    }
729    m_cache_test.erase(m_cache_test.begin(), m_cache_test.end());
730#endif   // CACHE
731
732    result->fullres = res;
733    result->length = length+1;
734    return result;
735}
736
Note: See TracBrowser for help on using the repository browser.