source: sage/libs/singular/singular.pyx @ 5614:cfee025945bc

Revision 5614:cfee025945bc, 7.7 KB checked in by 'Martin Albrecht <malb@…, 6 years ago (diff)

fixed SEGFAULT in libSINGULAR conversion routines.

Testcase:

sage: K = GF(next_prime(104)2, 'a')
sage: R.<x,y,z> = K[]
sage: z._rmul_(K(1))

Reported by Robert Bradshaw
Line 
1"""
2Singular C function and class declaration
3
4AUTHOR: Martin Albrecht <malb@informatik.uni-bremen.de>
5"""
6
7################################################################################
8#
9################################################################################
10
11###############################################################################
12#   SAGE: System for Algebra and Geometry Experimentation   
13#       Copyright (C) 2005, 2006 William Stein <wstein@gmail.com>
14#  Distributed under the terms of the GNU General Public License (GPL)
15#  The full text of the GPL is available at:
16#                  http://www.gnu.org/licenses/
17###############################################################################
18
19include "singular-cdefs.pxi"
20
21cdef extern from "limits.h":
22    long INT_MAX
23    long INT_MIN
24
25from sage.rings.rational_field import RationalField
26from sage.rings.finite_field import FiniteField_prime_modn
27from sage.rings.finite_field import FiniteField_ext_pari
28from sage.libs.pari.all import pari
29
30cdef extern from "stdsage.h":
31    ctypedef void PyObject
32    object PY_NEW(object t)
33    int PY_TYPE_CHECK(object o, object t)
34    PyObject** FAST_SEQ_UNSAFE(object o)
35    void init_csage()
36
37cdef class Conversion:
38    """
39    A work around class to export the contained methods/functions
40    """
41
42    cdef public Rational si2sa_QQ(self, number *n, ring *_ring):
43        """
44        Converts a SINGULAR rational number to a SAGE rational number.
45
46        INPUT:
47            n -- number
48            _ring -- singular ring, used to check type of n
49
50        OUTPUT:
51            SAGE rational number matching n
52        """
53
54        #TYPECHECK HERE
55
56        cdef number *nom
57        cdef number *denom
58        cdef mpq_t _z
59
60        cdef mpz_t nom_z, denom_z
61
62        cdef Rational z
63
64        mpq_init(_z)
65
66        ##  Immediate integers handles carry the tag 'SR_INT', i.e. the last bit is 1.
67        ##  This distuingishes immediate integers from other handles which  point  to
68        ##  structures aligned on 4 byte boundaries and therefor have last bit  zero.
69        ##  (The second bit is reserved as tag to allow extensions of  this  scheme.)
70        ##  Using immediates as pointers and dereferencing them gives address errors.
71       
72        nom = nlGetNom(n, _ring)
73        mpz_init(nom_z)
74
75        if (SR_HDL(nom) & SR_INT): mpz_set_si(nom_z, SR_TO_INT(nom))
76        else: mpz_set(nom_z,&nom.z)
77
78        mpq_set_num(_z,nom_z)
79        n_Delete(&nom,_ring)
80        mpz_clear(nom_z)
81
82        denom = nlGetDenom(n, _ring)
83        mpz_init(denom_z)
84
85        if (SR_HDL(denom) & SR_INT): mpz_set_si(denom_z, SR_TO_INT(denom))
86        else: mpz_set(denom_z,&denom.z)
87
88        mpq_set_den(_z, denom_z)
89        n_Delete(&denom,_ring)
90        mpz_clear(denom_z)
91
92        z = Rational()
93        z.set_from_mpq(_z)
94        return z
95
96    cdef public FiniteField_givaroElement si2sa_GFqGivaro(self, number *n, ring *_ring, FiniteField_givaro base):
97        cdef napoly *z
98        cdef int c, e
99        cdef int a
100        cdef int ret
101
102        if naIsZero(n):
103            return base._zero_element
104        elif naIsOne(n):
105            return base._one_element
106        z = (<lnumber*>n).z
107
108        a = base.objectptr.sage_generator()
109        ret = base.objectptr.zero
110
111        while z:
112            c = base.objectptr.read(c,<long>napGetCoeff(z))
113            e = napGetExp(z,1)
114            if e == 0:
115                ret = base.objectptr.add(ret, <int>c, ret)
116            else:
117                a = e * base.objectptr.sage_generator()
118                ret = base.objectptr.axpy(ret, <int>c, a, ret)
119            z = napIter(z)
120        return (<FiniteField_givaroElement>base._zero_element)._new_c(ret)
121           
122    cdef public object si2sa_GFqPari(self, number *n, ring *_ring, object base):
123        cdef napoly *z
124        cdef int c, e
125        cdef object a
126        cdef object ret
127
128        if naIsZero(n):
129            return base._zero_element
130        elif naIsOne(n):
131            return base._one_element
132        z = (<lnumber*>n).z
133
134        a = pari("a")
135        ret = pari(int(0)).Mod(int(_ring.ch))
136
137        while z:
138            c = <long>napGetCoeff(z)
139            e = napGetExp(z,1)
140            if e == 0:
141                ret = ret + c
142            elif c != 0:
143                ret = ret  + c * a**e
144            z = napIter(z)
145        return base(ret)
146
147
148
149    cdef public number *sa2si_QQ(self, Rational r, ring *_ring):
150        """
151        """
152        return nlInit2gmp( mpq_numref(r.value), mpq_denref(r.value) )
153
154    cdef number *sa2si_GFqGivaro(self, int quo, ring *_ring):
155        """
156        """
157        #can be done much faster
158        cdef number *n1, *n2, *a, *coeff, *apow1, *apow2
159        cdef int b
160       
161        rChangeCurrRing(_ring)
162        b   = - _ring.ch;
163       
164        a = naPar(1)
165
166        apow1 = naInit(1)
167        n1 = naInit(0)
168       
169        while quo!=0:
170            coeff = naInit(quo%b)
171           
172            if not naIsZero(coeff):
173                n2 = naAdd( naMult(coeff, apow1),  n1)
174                naDelete(&n1, _ring);
175                n1= n2
176
177            apow2 = naMult(apow1, a)
178            naDelete(&apow1, _ring)
179            apow1 = apow2
180           
181            quo = quo/b
182            naDelete(&coeff, _ring)
183           
184        naDelete(&apow1, _ring)
185        naDelete(&a, _ring)
186        return n1
187
188    cdef number *sa2si_GFqPari(self, object elem, ring *_ring):
189        #can be done much faster
190        cdef int i
191        cdef number *n1, *n2, *a, *coeff, *apow1, *apow2
192       
193        rChangeCurrRing(_ring)
194
195        elem = elem._pari_().lift().lift()
196       
197
198        if len(elem) > 1:
199            n1 = naInit(0)
200            a = naPar(1)
201            apow1 = naInit(1)
202
203            for i from 0 <= i < len(elem):
204                coeff = naInit(int(elem[i]))
205
206                if not naIsZero(coeff):
207                    n2 = naAdd( naMult(coeff, apow1),  n1)
208                    naDelete(&n1, _ring);
209                    n1= n2
210
211                apow2 = naMult(apow1, a)
212                naDelete(&apow1, _ring)
213                apow1 = apow2
214
215                naDelete(&coeff, _ring)
216
217            naDelete(&apow1, _ring)
218            naDelete(&a, _ring)
219        else:
220            n1 = naInit(int(elem))
221           
222        return n1
223
224    cdef public number *sa2si_ZZ(self, Integer d, ring *_ring):
225        """
226        """
227        cdef number *n
228        if INT_MIN <= d <= INT_MAX:
229            return nlInit(int(d))
230        else:
231            n = nlRInit(0)
232            mpz_init_set(&n.z, d.value)
233            return n
234
235    cdef public object si2sa(self, number *n, ring *_ring, object base):
236        if PY_TYPE_CHECK(base, FiniteField_prime_modn):
237            return base(nInt(n))
238           
239        elif PY_TYPE_CHECK(base, RationalField):
240            return self.si2sa_QQ(n,_ring)
241
242        elif PY_TYPE_CHECK(base, FiniteField_givaro):
243            return self.si2sa_GFqGivaro(n, _ring, base)
244
245        elif PY_TYPE_CHECK(base, FiniteField_ext_pari):
246            return self.si2sa_GFqPari(n, _ring, base)
247       
248        else:
249            raise ValueError, "cannot convert from SINGULAR number"
250
251    cdef public number *sa2si(self, Element elem, ring * _ring):
252        cdef int i
253        if PY_TYPE_CHECK(elem._parent, FiniteField_prime_modn):
254            return n_Init(int(elem),_ring)
255           
256        elif PY_TYPE_CHECK(elem._parent, RationalField):
257            return self.sa2si_QQ(elem, _ring)
258
259        elif PY_TYPE_CHECK(elem._parent, FiniteField_givaro):
260            return self.sa2si_GFqGivaro( (<FiniteField_givaro>elem._parent).objectptr.write(i, (<FiniteField_givaroElement>elem).element ), _ring )
261
262        elif PY_TYPE_CHECK(elem._parent, FiniteField_ext_pari):
263            return self.sa2si_GFqPari(elem, _ring)
264        else:
265            raise ValueError, "cannot convert to SINGULAR number"
266
Note: See TracBrowser for help on using the repository browser.