1 | from Singular import * |
---|
2 | from interpreter import * |
---|
3 | from objects import * |
---|
4 | from util import create_ring |
---|
5 | import CD.polyd |
---|
6 | from CD.polyd import DMPsym, SDMPsym, termsym, lpsym,dpsym,poly_ring_dsym, DMPLsym, implementation, groebnerdsym |
---|
7 | from CD.fieldname1 import Qsym as Rationals |
---|
8 | singular=singular_globals_proxy() |
---|
9 | class SingularException(Exception): |
---|
10 | pass |
---|
11 | |
---|
12 | def encodePoly(p): |
---|
13 | terms=[encodeTerm(t) for t in p] |
---|
14 | return OMApply(SDMPsym,terms) |
---|
15 | |
---|
16 | def encodePolyWithRing(p): |
---|
17 | pe=encodePoly(p) |
---|
18 | r=encodeRing(p.ring()) |
---|
19 | return OMApply(DMPsym,[r,pe]) |
---|
20 | |
---|
21 | def encodeGB(gb): |
---|
22 | i=encodeIdeal(gb) |
---|
23 | o=encodeOrdering(gb.ring()) |
---|
24 | return OMApply(groebnerdsym,[o,i]) |
---|
25 | orderingTable={ |
---|
26 | "lp": lpsym, |
---|
27 | "dp": dpsym |
---|
28 | } |
---|
29 | OrderingTableBack={} |
---|
30 | for k in orderingTable: |
---|
31 | OrderingTableBack[orderingTable[k]]=k |
---|
32 | |
---|
33 | def encodeOrdering(r): |
---|
34 | rl=singular.ringlist(r) |
---|
35 | return orderingTable[rl[2][0][0]] |
---|
36 | |
---|
37 | def encodeField(r): |
---|
38 | char=singular.char(r) |
---|
39 | if char==0 and singular.npars(r)==0: |
---|
40 | return Rationals |
---|
41 | else: |
---|
42 | raise SingularException("unknown field to encode") |
---|
43 | |
---|
44 | def encodeIdeal(i): |
---|
45 | r=encodeRing(i.ring()) |
---|
46 | return OMApply(DMPLsym,[r]+[encodePoly(p) for p in i]) |
---|
47 | def encodeRing(r): |
---|
48 | nv=singular.nvars(r) |
---|
49 | f=encodeField(r) |
---|
50 | return OMApply(poly_ring_dsym,[f,OMint(nv)]) |
---|
51 | |
---|
52 | def ringFromOM(ring_desc, ordering="dp"): |
---|
53 | assert isinstance(ring_desc, OMApply) |
---|
54 | if (ring_desc.args[0]==Rationals): |
---|
55 | i=ring_desc.args[1].getValue() |
---|
56 | return create_ring(char=0, nvars=i, ordering=ordering) |
---|
57 | raise SingularException("ring not supported") |
---|
58 | def polyFromOM(poly_desc): |
---|
59 | """assumes the right ring is set""" |
---|
60 | assert isinstance(poly_desc, OMApply) |
---|
61 | terms=[termFromOM(t) for t in poly_desc.args] |
---|
62 | res=polynomial(number(0)) |
---|
63 | for t in terms: |
---|
64 | res=res+t |
---|
65 | #res+=t |
---|
66 | return res |
---|
67 | def termFromOM(term_desc): |
---|
68 | assert isinstance(term_desc, OMApply) |
---|
69 | assert len(term_desc.args)==singular.nvars(ring())+1 |
---|
70 | assert isinstance(term_desc.args[0], OMint) |
---|
71 | |
---|
72 | coef=number(term_desc.args[0].getValue()) |
---|
73 | exp=intvec() |
---|
74 | for e in term_desc.args[1:]: |
---|
75 | assert isinstance(e,OMint) |
---|
76 | exp.append(e.getValue()) |
---|
77 | #print coef, polynomial(exp) |
---|
78 | return coef*polynomial(exp) |
---|
79 | |
---|
80 | def idealFromDMPL(dmpl): |
---|
81 | """assumes that the right ring is set""" |
---|
82 | i=ideal() |
---|
83 | assert len(dmpl.args)>=1 |
---|
84 | ps=[polyFromOM(d) for d in dmpl.args[1:]] |
---|
85 | for p in ps: |
---|
86 | i.append(p) |
---|
87 | return i |
---|
88 | def ringFromDMPLOrd(dmpl,o): |
---|
89 | assert len(dmpl.args)>=1 |
---|
90 | return ringFromOM(dmpl.args[0], ordering=OrderingTableBack[o]) |
---|
91 | leadcoef=singular.leadcoef |
---|
92 | leadexp=singular.leadexp |
---|
93 | def encodeTerm(t): |
---|
94 | """FIXME: ugly because it uses slow interpreter interface and setting of rings for this should be automatically""" |
---|
95 | #t.ring().set() |
---|
96 | #exponents=leadexp(t) |
---|
97 | #c=leadcoef(t) |
---|
98 | exponents=t.leadExp() |
---|
99 | c=t.leadCoef() |
---|
100 | exponents=[OMint(i) for i in exponents] |
---|
101 | return OMApply(termsym,[OMint(str(c))]+exponents) |
---|
102 | def groebnerfunc(context, ordering, dmpl): |
---|
103 | r=ringFromDMPLOrd(dmpl,ordering) |
---|
104 | r.set() |
---|
105 | i=idealFromDMPL(dmpl) |
---|
106 | res=singular.std(i) |
---|
107 | #FIXME: singular.groebner does not work may because of bad ring changes |
---|
108 | return encodeGB(res) |
---|
109 | implementation.implement("groebner", groebnerfunc) |
---|
110 | optimize(encodePoly) |
---|
111 | optimize(encodeTerm) |
---|