Context Navigation

source: sage/rings/polynomial/multi_polynomial_libsingular.pyx @ 5515:df32561a3872

Visit:

Revision 5515:df32561a3872, 106.4 KB checked in by 'Martin Albrecht <malb@…, 6 years ago (diff)
faster f.subs() for MPolynomial_libsingular

Line
1	"""
2	Multivariate polynomials via libSINGULAR.
3
4	AUTHORS:
5	-- Martin Albrecht <malb@informatik.uni-bremen.de>
6
7	TODO:
8	-- implement $GF(p^n)$
9	-- implement block orderings
10	-- implement Real, Complex
11
12	TESTS:
13	sage: from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
14	sage: P.<x,y,z> = MPolynomialRing_libsingular(QQ,3)
15	sage: loads(dumps(P)) == P
16	True
17	sage: loads(dumps(x)) == x
18	True
19	sage: P.<x,y,z> = MPolynomialRing_libsingular(GF(2^8,'a'),3)
20	sage: loads(dumps(P)) == P
21	True
22	sage: loads(dumps(x)) == x
23	True
24	sage: P.<x,y,z> = MPolynomialRing_libsingular(GF(127),3)
25	sage: loads(dumps(P)) == P
26	True
27	sage: loads(dumps(x)) == x
28	True
29
30	"""
31
32	include "sage/ext/interrupt.pxi"
33
34	cdef extern from "stdsage.h":
35	ctypedef void PyObject
36	object PY_NEW(object t)
37	int PY_TYPE_CHECK(object o, object t)
38	PyObject** FAST_SEQ_UNSAFE(object o)
39	void init_csage()
40
41	void sage_free(void *p)
42	void* sage_realloc(void *p, size_t n)
43	void* sage_malloc(size_t)
44
45	import os
46	import sage.rings.memory
47
48	from sage.libs.singular.singular import Conversion
49	from sage.libs.singular.singular cimport Conversion
50
51	cdef Conversion co
52	co = Conversion()
53
54	from sage.rings.polynomial.multi_polynomial_ring import MPolynomialRing_polydict_domain
55	from sage.rings.polynomial.term_order import TermOrder
56	from sage.rings.polynomial.multi_polynomial_element import MPolynomial_polydict
57	from sage.rings.polynomial.multi_polynomial_ideal import MPolynomialIdeal
58	from sage.rings.polynomial.polydict import ETuple
59	from sage.rings.polynomial.polynomial_ring import PolynomialRing
60
61	from sage.rings.rational_field import RationalField
62	from sage.rings.finite_field import FiniteField_prime_modn
63	from sage.rings.finite_field import FiniteField_generic
64	from sage.rings.finite_field_givaro cimport FiniteField_givaro
65
66	from sage.rings.number_field.number_field import NumberField_generic
67
68	from sage.rings.finite_field_givaro cimport FiniteField_givaroElement
69
70	from sage.rings.rational cimport Rational
71
72	from sage.interfaces.singular import singular as singular_default, is_SingularElement, SingularElement
73	from sage.interfaces.macaulay2 import macaulay2 as macaulay2_default, is_Macaulay2Element
74	from sage.structure.factorization import Factorization
75
76	from sage.rings.complex_field import is_ComplexField
77	from sage.rings.real_mpfr import is_RealField
78	from sage.rings.integer_ring import is_IntegerRing
79
80	from sage.rings.integer_ring import IntegerRing
81	from sage.structure.element cimport EuclideanDomainElement, \
82	RingElement, \
83	ModuleElement, \
84	Element, \
85	CommutativeRingElement
86
87	from sage.structure.sequence import Sequence
88
89
90	from sage.rings.integer cimport Integer
91
92	from sage.structure.parent cimport Parent
93	from sage.structure.parent_base cimport ParentWithBase
94	from sage.structure.parent_gens cimport ParentWithGens
95
96	from sage.misc.misc import mul
97	from sage.misc.sage_eval import sage_eval
98	from sage.misc.latex import latex
99
100	from sage.interfaces.all import macaulay2
101
102
103	# shared library loading
104	cdef extern from "dlfcn.h":
105	void dlopen(char , long)
106	char *dlerror()
107
108	cdef extern from "string.h":
109	char strdup(char s)
110
111	cdef init_singular():
112	"""
113	This initializes the Singular library. Right now, this is a hack.
114
115	SINGULAR has a concept of compiled extension modules similar to
116	SAGE. For this, the compiled modules need to see the symbols from
117	the main programm. However, SINGULAR is a shared library in this
118	context these symbols are not known globally. The work around so
119	far is to load the library again and to specifiy RTLD_GLOBAL.
120	"""
121	cdef void *handle
122
123
124	for extension in ["so", "dylib", "dll"]:
125	lib = os.environ['SAGE_LOCAL']+"/lib/libsingular."+extension
126	if os.path.exists(lib):
127	handle = dlopen(lib, 256+1)
128	break
129
130	if handle == NULL:
131	print dlerror()
132	raise ImportError, "cannot load libSINGULAR library"
133
134	# Load Singular
135	siInit(lib)
136
137	# Steal Memory Manager back or weird things may happen
138	sage.rings.memory.pmem_malloc()
139
140	# call it
141	init_singular()
142
143	order_dict = {"dp":ringorder_dp,
144	"Dp":ringorder_Dp,
145	"lp":ringorder_lp,
146	"rp":ringorder_rp,
147	"ds":ringorder_ds,
148	"Ds":ringorder_Ds,
149	"ls":ringorder_ls,
150	}
151
152	cdef class MPolynomialRing_libsingular(MPolynomialRing_generic):
153	"""
154	A multivariate polynomial ring over QQ or GF(p) implemented using SINGULAR.
155
156	"""
157	def __init__(self, base_ring, n, names, order='degrevlex'):
158	"""
159
160	Construct a multivariate polynomial ring subject to the following conditions:
161
162	INPUT:
163	base_ring -- base ring (must be either GF(p) (p prime) or QQ)
164	n -- number of variables (must be at least 1)
165	names -- names of ring variables, may be string of list/tuple
166	order -- term order (default: degrevlex)
167
168	EXAMPLES:
169	sage: P.<x,y,z> = QQ[]
170	sage: P
171	Polynomial Ring in x, y, z over Rational Field
172
173	sage: f = 27/113 * x^2 + y*z + 1/2; f
174	27/113x^2 + yz + 1/2
175
176	sage: P.term_order()
177	Degree reverse lexicographic term order
178
179	sage: P = MPolynomialRing(GF(127),3,names='abc', order='lex')
180	sage: P
181	Polynomial Ring in a, b, c over Finite Field of size 127
182
183	sage: a,b,c = P.gens()
184	sage: f = 57 * a^2b + 43 c + 1; f
185	57a^2b + 43*c + 1
186
187	sage: P.term_order()
188	Lexicographic term order
189
190	"""
191	cdef char **_names
192	cdef char *_name
193	cdef int i
194	cdef int nblcks
195	cdef int offset
196	cdef int characteristic
197	cdef MPolynomialRing_libsingular k
198	cdef MPolynomial_libsingular minpoly
199	cdef lnumber *nmp
200
201	is_extension = False
202
203	n = int(n)
204	if n<1:
205	raise ArithmeticError, "number of variables must be at least 1"
206
207	self.__ngens = n
208
209	order = TermOrder(order, n)
210
211	MPolynomialRing_generic.__init__(self, base_ring, n, names, order)
212
213	self._has_singular = True
214
215	_names = <char*>omAlloc0(sizeof(char)*(len(self._names)+1))
216
217	for i from 0 <= i < n:
218	_name = self._names[i]
219	_names[i] = omStrDup(_name)
220
221	# from the SINGULAR source code documentation for the rInit function
222	## characteristic --------------------------------------------------
223	## input: 0 ch=0 : Q parameter=NULL ffChar=FALSE float_len (done)
224	## 0 1 : Q(a,...) *names FALSE (todo)
225	## 0 -1 : R NULL FALSE 0
226	## 0 -1 : R NULL FALSE prec. >6
227	## 0 -1 : C *names FALSE prec. 0..?
228	## p p : Fp NULL FALSE (done)
229	## p -p : Fp(a) *names FALSE (done)
230	## q q : GF(q=p^n) *names TRUE (todo)
231
232	if PY_TYPE_CHECK(base_ring, FiniteField_prime_modn):
233	if base_ring.characteristic() <= 2147483629:
234	characteristic = base_ring.characteristic()
235	else:
236	raise TypeError, "p must be <= 2147483629"
237
238	elif PY_TYPE_CHECK(base_ring, RationalField):
239	characteristic = 0
240
241	elif PY_TYPE_CHECK(base_ring, FiniteField_generic):
242	characteristic = -base_ring.characteristic() # note the negative characteristic
243	k = MPolynomialRing_libsingular(base_ring.prime_subfield(), 1, base_ring.variable_name(), 'lex')
244	minpoly = base_ring.polynomial()(k.gen())
245	is_extension = True
246
247	elif PY_TYPE_CHECK(base_ring, NumberField_generic):
248	characteristic = 1
249	k = MPolynomialRing_libsingular(RationalField(), 1, base_ring.variable_name(), 'lex')
250	minpoly = base_ring.polynomial()(k.gen())
251	is_extension = True
252
253	else:
254	raise NotImplementedError, "Only GF(q), QQ, and Number Fields are supported."
255
256	self._ring = <ring*>omAlloc0Bin(sip_sring_bin)
257	self._ring.ch = characteristic
258	self._ring.N = n
259	self._ring.names = _names
260
261	if is_extension:
262	rChangeCurrRing(k._ring)
263	self._ring.algring = rCopy0(k._ring)
264	rComplete(self._ring.algring,1)
265	self._ring.P = self._ring.algring.N
266	#self._ring.parameter = self._ring.algring.names
267	self._ring.parameter = <char*>omAlloc0(sizeof(char)*2)
268	self._ring.parameter[0] = omStrDup(self._ring.algring.names[0])
269
270	nmp = <lnumber*>omAlloc0Bin(rnumber_bin)
271	nmp.z= <napoly*>p_Copy(minpoly._poly, self._ring.algring) # fragile?
272	nmp.s=2
273
274	self._ring.minpoly=<number*>nmp
275
276	nblcks = len(order.blocks)
277	offset = 0
278
279	self._ring.wvhdl = <int *>omAlloc0((nblcks + 2) sizeof(int*))
280	self._ring.order = <int >omAlloc0((nblcks + 2) sizeof(int *))
281	self._ring.block0 = <int >omAlloc0((nblcks + 2) sizeof(int *))
282	self._ring.block1 = <int >omAlloc0((nblcks + 2) sizeof(int *))
283	self._ring.OrdSgn = 1
284
285
286	for i from 0 <= i < nblcks:
287	self._ring.order[i] = order_dict.get(order.blocks[i][0], ringorder_lp)
288	self._ring.block0[i] = offset + 1
289	if order.blocks[i][1] == 0: # may be zero in some cases
290	self._ring.block1[i] = offset + n
291	else:
292	self._ring.block1[i] = offset + order.blocks[i][1]
293	offset = self._ring.block1[i]
294
295	self._ring.order[nblcks] = ringorder_C
296
297	rComplete(self._ring, 1)
298	self._ring.ShortOut = 0
299
300	self._zero = <MPolynomial_libsingular>new_MP(self,NULL)
301
302	def __dealloc__(self):
303	"""
304	"""
305	rChangeCurrRing(self._ring)
306	rDelete(self._ring)
307
308	cdef _coerce_c_impl(self, element):
309	"""
310
311	Coerces elements to self.
312
313	EXAMPLES:
314	sage: P.<x,y,z> = QQ[]
315
316	We can coerce elements of self to self
317
318	sage: P._coerce_(x*y + 1/2)
319	x*y + 1/2
320
321	We can coerce elements for a ring with the same algebraic properties
322
323	sage: from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
324	sage: R.<x,y,z> = MPolynomialRing_libsingular(QQ,3)
325	sage: P == R
326	True
327
328	sage: P is R
329	False
330
331	sage: P._coerce_(x*y + 1)
332	x*y + 1
333
334	We can coerce base ring elements
335
336	sage: P._coerce_(3/2)
337	3/2
338
339	sage: P._coerce_(ZZ(1))
340	1
341
342	sage: P._coerce_(int(1))
343	1
344
345	sage: k.<a> = GF(2^8)
346	sage: P.<x,y> = PolynomialRing(k,2)
347	sage: P._coerce_(a)
348	(a)
349
350	"""
351	cdef poly *_p
352	cdef ring *_ring
353	cdef number *_n
354	cdef poly *mon
355	cdef int i
356
357	_ring = self._ring
358
359	if(_ring != currRing): rChangeCurrRing(_ring)
360
361	if PY_TYPE_CHECK(element, MPolynomial_libsingular):
362	if element.parent() is <object>self:
363	return element
364	elif element.parent() == self:
365	# is this safe?
366	_p = p_Copy((<MPolynomial_libsingular>element)._poly, _ring)
367	else:
368	raise TypeError, "parents do not match"
369
370	elif PY_TYPE_CHECK(element, MPolynomial_polydict):
371	if element.parent() == self:
372	_p = p_ISet(0, _ring)
373	for (m,c) in element.element().dict().iteritems():
374	mon = p_Init(_ring)
375	p_SetCoeff(mon, co.sa2si(c, _ring), _ring)
376	for pos in m.nonzero_positions():
377	p_SetExp(mon, pos+1, m[pos], _ring)
378	p_Setm(mon, _ring)
379	_p = p_Add_q(_p, mon, _ring)
380	else:
381	raise TypeError, "parents do not match"
382
383	elif PY_TYPE_CHECK(element, CommutativeRingElement):
384	# base ring elements
385	if <Parent>element.parent() is self._base:
386	# shortcut for GF(p)
387	if PY_TYPE_CHECK(self._base, FiniteField_prime_modn):
388	_p = p_ISet(int(element), _ring)
389	else:
390	_n = co.sa2si(element,_ring)
391	_p = p_NSet(_n, _ring)
392
393	# also accepting ZZ
394	elif element.parent() is IntegerRing():
395	if PY_TYPE_CHECK(self._base, RationalField):
396	_n = co.sa2si_ZZ(element,_ring)
397	_p = p_NSet(_n, _ring)
398	else: # GF(p)
399	_p = p_ISet(int(element),_ring)
400	else:
401	# fall back to base ring
402	return self._base._coerce_c(element)
403
404	# Accepting int
405	elif PY_TYPE_CHECK(element, int):
406	_p = p_ISet(int(element), _ring)
407	else:
408	raise TypeError, "Cannot coerce element"
409
410	return new_MP(self,_p)
411
412	def __call__(self, element):
413	"""
414	Construct a new element in self.
415
416	INPUT:
417	element -- several types are supported, see below
418
419	EXAMPLE:
420	Call supports all conversions _coerce_ supports, plus:
421
422	Coercion form strings:
423	sage: P.<x,y,z> = QQ[]
424	sage: P('x+y + 1/4')
425	x + y + 1/4
426
427	Coercion from SINGULAR elements:
428	sage: P._singular_()
429	// characteristic : 0
430	// number of vars : 3
431	// block 1 : ordering dp
432	// : names x y z
433	// block 2 : ordering C
434
435	sage: P._singular_().set_ring()
436	sage: P(singular('x + 3/4'))
437	x + 3/4
438
439	Coercion from symbolic variables:
440	sage: x,y,z = var('x,y,z')
441	sage: R = QQ[x,y,z]
442	sage: R(x)
443	x
444
445	Coercion from 'similar' rings:
446	sage: P.<x,y,z> = QQ[]
447	sage: R.<a,b,c> = MPolynomialRing(ZZ,3)
448	sage: P(a)
449	x
450
451	If everything else fails, we try to coerce to the base ring:
452	sage: R.<x,y,z> = GF(3)[]
453	sage: R(1/2)
454	-1
455
456	"""
457	cdef poly _p, mon
458	cdef ring *_ring = self._ring
459	rChangeCurrRing(_ring)
460
461	# try to coerce first
462	try:
463	return self._coerce_c_impl(element)
464	except TypeError:
465	pass
466
467	if PY_TYPE_CHECK(element, SingularElement):
468	element = str(element)
469
470	if PY_TYPE_CHECK(element, basestring):
471	# let python do the the parsing
472	d = self.gens_dict()
473	if PY_TYPE_CHECK(self._base, FiniteField_givaro):
474	d[str(self._base.gen())]=self._base.gen()
475	element = sage_eval(element,d)
476
477	# we need to do this, to make sure that we actually get an
478	# element in self.
479	return self._coerce_c(element)
480
481	if PY_TYPE_CHECK(element, MPolynomial_libsingular):
482	if element.parent() is not self and element.parent() != self and element.parent().ngens() == self.ngens():
483	# Map the variables in some crazy way (but in order,
484	# of course). This is here since R(blah) is supposed
485	# to be "make an element of R if at all possible with
486	# no guarantees that this is mathematically solid."
487	# TODO: We can do this faster without the dict
488	_p = p_ISet(0, _ring)
489	K = self.base_ring()
490	for (m,c) in element.dict().iteritems():
491	try:
492	c = K(c)
493	except TypeError, msg:
494	p_Delete(&_p, _ring)
495	raise TypeError, msg
496	mon = p_Init(_ring)
497	p_SetCoeff(mon, co.sa2si(c , _ring), _ring)
498	for pos in m.nonzero_positions():
499	p_SetExp(mon, pos+1, m[pos], _ring)
500	p_Setm(mon, _ring)
501	_p = p_Add_q(_p, mon, _ring)
502	return new_MP(self, _p)
503
504	if PY_TYPE_CHECK(element, MPolynomial_polydict):
505	if element.parent().ngens() == self.ngens():
506	# Map the variables in some crazy way (but in order,
507	# of course). This is here since R(blah) is supposed
508	# to be "make an element of R if at all possible with
509	# no guarantees that this is mathematically solid."
510	_p = p_ISet(0, _ring)
511	K = self.base_ring()
512	for (m,c) in element.element().dict().iteritems():
513	try:
514	c = K(c)
515	except TypeError, msg:
516	p_Delete(&_p, _ring)
517	raise TypeError, msg
518	mon = p_Init(_ring)
519	p_SetCoeff(mon, co.sa2si(c , _ring), _ring)
520	for pos in m.nonzero_positions():
521	p_SetExp(mon, pos+1, m[pos], _ring)
522	p_Setm(mon, _ring)
523	_p = p_Add_q(_p, mon, _ring)
524	return new_MP(self, _p)
525
526	if hasattr(element,'_polynomial_'):
527	# SymbolicVariable
528	return element._polynomial_(self)
529
530	if is_Macaulay2Element(element):
531	return self(repr(element))
532
533	# now try calling the base ring's __call__ methods
534	element = self.base_ring()(element)
535	_p = p_NSet(co.sa2si(element,_ring), _ring)
536	return new_MP(self,_p)
537
538	def _repr_(self):
539	"""
540	EXAMPLE:
541	sage: P.<x,y> = QQ[]
542	sage: P
543	Polynomial Ring in x, y over Rational Field
544
545	"""
546	varstr = ", ".join([ rRingVar(i,self._ring) for i in range(self.__ngens) ])
547	return "Polynomial Ring in %s over %s"%(varstr,self._base)
548
549	def ngens(self):
550	"""
551	Returns the number of variables in self.
552
553	EXAMPLES:
554	sage: P.<x,y> = QQ[]
555	sage: P.ngens()
556	2
557
558	sage: k.<a> = GF(2^16)
559	sage: P = PolynomialRing(k,1000,'x')
560	sage: P.ngens()
561	1000
562
563	"""
564	return int(self.__ngens)
565
566	def gens(self):
567	"""
568	Return the tuple of variables in self.
569
570	EXAMPLES:
571	sage: P.<x,y,z> = QQ[]
572	sage: P.gens()
573	(x, y, z)
574
575	sage: P = MPolynomialRing(QQ,10,'x')
576	sage: P.gens()
577	(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9)
578
579	sage: P.<SAGE,SINGULAR> = MPolynomialRing(QQ,2) # weird names
580	sage: P.gens()
581	(SAGE, SINGULAR)
582
583	"""
584	return tuple([self.gen(i) for i in range(self.__ngens) ])
585
586	def gen(self, int n=0):
587	"""
588	Returns the n-th generator of self.
589
590	EXAMPLES:
591	sage: P.<x,y,z> = QQ[]
592	sage: P.gen(),P.gen(1)
593	(x, y)
594
595	sage: P = MPolynomialRing(GF(127),1000,'x')
596	sage: P.gen(500)
597	x500
598
599	sage: P.<SAGE,SINGULAR> = MPolynomialRing(QQ,2) # weird names
600	sage: P.gen(1)
601	SINGULAR
602
603	"""
604	cdef poly *_p
605	cdef ring *_ring = self._ring
606
607	if n < 0 or n >= self.__ngens:
608	raise ValueError, "Generator not defined."
609
610	rChangeCurrRing(_ring)
611	_p = p_ISet(1,_ring)
612	p_SetExp(_p, n+1, 1, _ring)
613	p_Setm(_p, _ring);
614
615	return new_MP(self,_p)
616
617	def ideal(self, gens, coerce=True):
618	"""
619	Create an ideal in this polynomial ring.
620
621	INPUT:
622	gens -- generators of the ideal
623	coerce -- shall the generators be coerced first (default:True)
624
625	EXAMPLE:
626	sage: P.<x,y,z> = QQ[]
627	sage: sage.rings.ideal.Katsura(P)
628	Ideal (x + 2y + 2z - 1, x^2 + 2y^2 + 2z^2 - x, 2xy + 2yz - y) of Polynomial Ring in x, y, z over Rational Field
629
630	sage: P.ideal([x + 2y + 2z-1, 2xy + 2yz-y, x^2 + 2y^2 + 2z^2-x])
631	Ideal (x + 2y + 2z - 1, 2xy + 2yz - y, x^2 + 2y^2 + 2z^2 - x) of Polynomial Ring in x, y, z over Rational Field
632
633	"""
634	if is_SingularElement(gens):
635	gens = list(gens)
636	coerce = True
637	if is_Macaulay2Element(gens):
638	gens = list(gens)
639	coerce = True
640	elif not isinstance(gens, (list, tuple)):
641	gens = [gens]
642	if coerce:
643	gens = [self(x) for x in gens] # this will even coerce from singular ideals correctly!
644	return MPolynomialIdeal(self, gens, coerce=False)
645
646	def _macaulay2_(self, macaulay2=macaulay2_default):
647	"""
648	Create a M2 representation of self if Macaulay2 is installed.
649
650	INPUT:
651	macaulay2 -- M2 interpreter (default: macaulay2_default)
652
653	EXAMPLES:
654	sage: R.<x,y> = ZZ[]
655	sage: macaulay2(R) # optional
656	ZZ [x, y, MonomialOrder => GRevLex, MonomialSize => 16]
657
658	sage: R.<x,y> = QQ[]
659	sage: macaulay2(R) # optional
660	QQ [x, y, MonomialOrder => GRevLex, MonomialSize => 16]
661
662	sage: R.<x,y> = GF(17)[]
663	sage: print macaulay2(R) # optional
664	ZZ
665	-- [x, y, MonomialOrder => GRevLex, MonomialSize => 16]
666	17
667	"""
668	try:
669	R = self.__macaulay2
670	if R is None or not (R.parent() is macaulay2):
671	raise ValueError
672	R._check_valid()
673	return R
674	except (AttributeError, ValueError):
675	self.__macaulay2 = self._macaulay2_set_ring(macaulay2)
676	return self.__macaulay2
677
678	def _macaulay2_set_ring(self, macaulay2):
679	if not self.__m2_set_ring_cache is None:
680	base_str, gens, order = self.__m2_set_ring_cache
681	else:
682	if self.base_ring().is_prime_field():
683	if self.characteristic() == 0:
684	base_str = "QQ"
685	else:
686	base_str = "ZZ/" + str(self.characteristic())
687	elif is_IntegerRing(self.base_ring()):
688	base_str = "ZZ"
689	else:
690	raise TypeError, "no conversion of to a Macaulay2 ring defined"
691	gens = str(self.gens())
692	order = self.term_order().macaulay2_str()
693	self.__m2_set_ring_cache = (base_str, gens, order)
694	return macaulay2.ring(base_str, gens, order)
695
696	def _can_convert_to_singular(self):
697	"""
698	Returns True
699	"""
700	return True
701
702	def _singular_(self, singular=singular_default):
703	"""
704	Create a SINGULAR (as in the CAS) representation of self. The
705	result is cached.
706
707	INPUT:
708	singular -- SINGULAR interpreter (default: singular_default)
709
710	EXAMPLES:
711	sage: P.<x,y,z> = QQ[]
712	sage: P._singular_()
713	// characteristic : 0
714	// number of vars : 3
715	// block 1 : ordering dp
716	// : names x y z
717	// block 2 : ordering C
718
719	sage: P._singular_() is P._singular_()
720	True
721
722	sage: P._singular_().name() == P._singular_().name()
723	True
724
725	sage: k.<a> = GF(3^3)
726	sage: P.<x,y,z> = PolynomialRing(k,3)
727	sage: P._singular_()
728	// characteristic : 3
729	// 1 parameter : a
730	// minpoly : (a^3-a+1)
731	// number of vars : 3
732	// block 1 : ordering dp
733	// : names x y z
734	// block 2 : ordering C
735
736	sage: P._singular_() is P._singular_()
737	True
738
739	sage: P._singular_().name() == P._singular_().name()
740	True
741
742
743	TESTS:
744	sage: from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
745	sage: P.<x> = MPolynomialRing_libsingular(QQ,1)
746	sage: P._singular_()
747	// characteristic : 0
748	// number of vars : 1
749	// block 1 : ordering lp
750	// : names x
751	// block 2 : ordering C
752
753	"""
754	try:
755	R = self.__singular
756	if R is None or not (R.parent() is singular):
757	raise ValueError
758	R._check_valid()
759	if self.base_ring().is_prime_field():
760	return R
761	if self.base_ring().is_finite():
762	R.set_ring() #sorry for that, but needed for minpoly
763	if singular.eval('minpoly') != self.__minpoly:
764	singular.eval("minpoly=%s"%(self.__minpoly))
765	return R
766	except (AttributeError, ValueError):
767	return self._singular_init_(singular)
768
769	def _singular_init_(self, singular=singular_default):
770	"""
771	Create a SINGULAR (as in the CAS) representation of self. The
772	result is NOT cached.
773
774	INPUT:
775	singular -- SINGULAR interpreter (default: singular_default)
776
777	EXAMPLES:
778	sage: P.<x,y,z> = QQ[]
779	sage: P._singular_init_()
780	// characteristic : 0
781	// number of vars : 3
782	// block 1 : ordering dp
783	// : names x y z
784	// block 2 : ordering C
785	sage: P._singular_init_() is P._singular_init_()
786	False
787
788	sage: P._singular_init_().name() == P._singular_init_().name()
789	False
790
791	TESTS:
792	sage: P.<x> = QQ[]
793	sage: P._singular_init_()
794	// characteristic : 0
795	// number of vars : 1
796	// block 1 : ordering lp
797	// : names x
798	// block 2 : ordering C
799
800	"""
801	if self.ngens()==1:
802	_vars = str(self.gen())
803	if "" in _vars: # 1.000...000x
804	_vars = _vars.split("*")[1]
805	order = 'lp'
806	else:
807	_vars = str(self.gens())
808	order = self.term_order().singular_str()
809
810	if is_RealField(self.base_ring()):
811	# singular converts to bits from base_10 in mpr_complex.cc by:
812	# size_t bits = 1 + (size_t) ((float)digits * 3.5);
813	precision = self.base_ring().precision()
814	digits = sage.rings.arith.ceil((2*precision - 2)/7.0)
815	self.__singular = singular.ring("(real,%d,0)"%digits, _vars, order=order)
816
817	elif is_ComplexField(self.base_ring()):
818	# singular converts to bits from base_10 in mpr_complex.cc by:
819	# size_t bits = 1 + (size_t) ((float)digits * 3.5);
820	precision = self.base_ring().precision()
821	digits = sage.rings.arith.ceil((2*precision - 2)/7.0)
822	self.__singular = singular.ring("(complex,%d,0,I)"%digits, _vars, order=order)
823
824	elif self.base_ring().is_prime_field():
825	self.__singular = singular.ring(self.characteristic(), _vars, order=order)
826
827	elif self.base_ring().is_finite(): #must be extension field
828	gen = str(self.base_ring().gen())
829	r = singular.ring( "(%s,%s)"%(self.characteristic(),gen), _vars, order=order)
830	self.__minpoly = "("+(str(self.base_ring().modulus()).replace("x",gen)).replace(" ","")+")"
831	singular.eval("minpoly=%s"%(self.__minpoly) )
832
833	self.__singular = r
834	else:
835	raise TypeError, "no conversion to a Singular ring defined"
836
837	return self.__singular
838
839	def __hash__(self):
840	"""
841	Return a hash for self, that is, a hash of the string representation of self
842
843	EXAMPLE:
844	sage: P.<x,y,z> = QQ[]
845	sage: hash(P)
846	-6257278808099690586 # 64-bit
847	-1767675994 # 32-bit
848	"""
849	return hash(self.__repr__())
850
851	def __richcmp__(left, right, int op):
852	return (<Parent>left)._richcmp(right, op)
853
854	cdef int _cmp_c_impl(left, Parent right) except -2:
855	"""
856	Multivariate polynomial rings are said to be equal if:
857	* their base rings match
858	* their generator names match
859	* their term orderings match
860
861	EXAMPLES:
862	sage: P.<x,y,z> = QQ[]
863	sage: R.<x,y,z> = QQ[]
864	sage: P == R
865	True
866
867	sage: R.<x,y,z> = MPolynomialRing(GF(127),3)
868	sage: P == R
869	False
870
871	sage: R.<x,y> = MPolynomialRing(QQ,2)
872	sage: P == R
873	False
874
875	sage: R.<x,y,z> = MPolynomialRing(QQ,3,order='revlex')
876	sage: P == R
877	False
878
879
880	"""
881	if PY_TYPE_CHECK(right, MPolynomialRing_libsingular) or PY_TYPE_CHECK(right, MPolynomialRing_polydict_domain):
882	return cmp( (left.base_ring(), map(str, left.gens()), left.term_order()),
883	(right.base_ring(), map(str, right.gens()), right.term_order())
884	)
885	else:
886	return cmp(type(left),type(right))
887
888	def __reduce__(self):
889	"""
890	Serializes self.
891
892	EXAMPLES:
893	sage: P.<x,y,z> = MPolynomialRing(QQ,3, order='degrevlex')
894	sage: P == loads(dumps(P))
895	True
896
897	sage: P = MPolynomialRing(GF(127),3,names='abc')
898	sage: P == loads(dumps(P))
899	True
900
901	sage: P = PolynomialRing(GF(2^8,'F'),3,names='abc')
902	sage: P == loads(dumps(P))
903	True
904
905	sage: P = PolynomialRing(GF(2^16,'B'),3,names='abc')
906	sage: P == loads(dumps(P))
907	True
908
909	"""
910	return sage.rings.polynomial.multi_polynomial_libsingular.unpickle_MPolynomialRing_libsingular, ( self.base_ring(),
911	map(str, self.gens()),
912	self.term_order() )
913
914
915	def __temporarily_change_names(self, names, latex_names):
916	"""
917	This is used by the variable names context manager.
918	"""
919	cdef ring *_ring = (<MPolynomialRing_libsingular>self)._ring
920	cdef char _names, _orig_names
921	cdef char *_name
922	cdef int i
923
924	if len(names) != _ring.N:
925	raise TypeError, "len(names) doesn't equal self.ngens()"
926
927	old = self._names, self._latex_names
928	(self._names, self._latex_names) = names, latex_names
929
930	_names = <char*>omAlloc0(sizeof(char)*_ring.N)
931	for i from 0 <= i < _ring.N:
932	_name = names[i]
933	_names[i] = omStrDup(_name)
934
935	_orig_names = _ring.names
936	_ring.names = _names
937
938	for i from 0 <= i < _ring.N:
939	omFree(_orig_names[i])
940	omFree(_orig_names)
941
942	return old
943
944	### The following methods are handy for implementing Groebner
945	### basis algorithms. They do only superficial type/sanity checks
946	### and should be called carefully.
947
948	def monomial_quotient(self, MPolynomial_libsingular f, MPolynomial_libsingular g, coeff=False):
949	"""
950	Return f/g, where both f and g are treated as
951	monomials. Coefficients are ignored by default.
952
953	INPUT:
954	f -- monomial
955	g -- monomial
956	coeff -- divide coefficents as well (default: False)
957
958	EXAMPLE:
959	sage: P.<x,y,z> = QQ[]
960	sage: P.monomial_quotient(3/2xy,x)
961	y
962
963	sage: P.monomial_quotient(3/2xy,x,coeff=True)
964	3/2*y
965
966	TESTS:
967	sage: R.<x,y,z> = QQ[]
968	sage: P.<x,y,z> = QQ[]
969	sage: P.monomial_quotient(x*y,x)
970	y
971
972	sage: P.monomial_quotient(x*y,R.gen())
973	y
974
975	sage: P.monomial_quotient(P(0),P(1))
976	0
977
978	sage: P.monomial_quotient(P(1),P(0))
979	Traceback (most recent call last):
980	...
981	ZeroDivisionError
982
983	sage: P.monomial_quotient(P(3/2),P(2/3), coeff=True)
984	9/4
985
986	sage: P.monomial_quotient(x,y) # Note the wrong result
987	xy^1048575z^1048575 # 64-bit
988	xy^65535z^65535 # 32-bit
989
990	sage: P.monomial_quotient(x,P(1))
991	x
992
993	NOTE: Assumes that the head term of f is a multiple of the
994	head term of g and return the multiplicant m. If this rule is
995	violated, funny things may happen.
996
997	"""
998	cdef poly *res
999	cdef ring *r = self._ring
1000
1001	if not <ParentWithBase>self is f._parent:
1002	f = self._coerce_c(f)
1003	if not <ParentWithBase>self is g._parent:
1004	g = self._coerce_c(g)
1005
1006	if(r != currRing): rChangeCurrRing(r)
1007
1008	if not f._poly:
1009	return self._zero
1010	if not g._poly:
1011	raise ZeroDivisionError
1012
1013	res = pDivide(f._poly,g._poly)
1014	if coeff:
1015	p_SetCoeff(res, n_Div( p_GetCoeff(f._poly, r) , p_GetCoeff(g._poly, r), r), r)
1016	else:
1017	p_SetCoeff(res, n_Init(1, r), r)
1018	return new_MP(self, res)
1019
1020	def monomial_is_divisible_by(self, MPolynomial_libsingular a, MPolynomial_libsingular b):
1021	"""
1022	Return False if b does not divide a and True otherwise.
1023
1024	INPUT:
1025	a -- monomial
1026	b -- monomial
1027
1028	EXAMPLES:
1029	sage: P.<x,y,z> = QQ[]
1030	sage: P.monomial_is_divisible_by(x^3y^2z^4, xyz)
1031	True
1032	sage: P.monomial_is_divisible_by(xyz, x^3y^2z^4)
1033	False
1034
1035	TESTS:
1036	sage: P.<x,y,z> = QQ[]
1037	sage: P.monomial_is_divisible_by(P(0),P(1))
1038	True
1039	sage: P.monomial_is_divisible_by(x,P(1))
1040	True
1041
1042	"""
1043	cdef poly *_a
1044	cdef poly *_b
1045	cdef ring *_r
1046	if a._parent is not b._parent:
1047	b = (<MPolynomialRing_libsingular>a._parent)._coerce_c(b)
1048
1049	_a = a._poly
1050	_b = b._poly
1051	_r = (<MPolynomialRing_libsingular>a._parent)._ring
1052
1053	if _b == NULL:
1054	raise ZeroDivisionError
1055	if _a == NULL:
1056	return True
1057
1058	if not p_DivisibleBy(_b, _a, _r):
1059	return False
1060	else:
1061	return True
1062
1063
1064	def monomial_lcm(self, MPolynomial_libsingular f, MPolynomial_libsingular g):
1065	"""
1066	LCM for monomials. Coefficients are ignored.
1067
1068	INPUT:
1069	f -- monomial
1070	g -- monomial
1071
1072	EXAMPLE:
1073	sage: P.<x,y,z> = QQ[]
1074	sage: P.monomial_lcm(3/2xy,x)
1075	x*y
1076
1077	TESTS:
1078	sage: R.<x,y,z> = QQ[]
1079	sage: P.<x,y,z> = QQ[]
1080	sage: P.monomial_lcm(x*y,R.gen())
1081	x*y
1082
1083	sage: P.monomial_lcm(P(3/2),P(2/3))
1084	1
1085
1086	sage: P.monomial_lcm(x,P(1))
1087	x
1088
1089	"""
1090	cdef poly *m = p_ISet(1,self._ring)
1091
1092	if not <ParentWithBase>self is f._parent:
1093	f = self._coerce_c(f)
1094	if not <ParentWithBase>self is g._parent:
1095	g = self._coerce_c(g)
1096
1097	if f._poly == NULL:
1098	if g._poly == NULL:
1099	return self._zero
1100	else:
1101	raise ArithmeticError, "cannot compute lcm of zero and nonzero element"
1102	if g._poly == NULL:
1103	raise ArithmeticError, "cannot compute lcm of zero and nonzero element"
1104
1105	if(self._ring != currRing): rChangeCurrRing(self._ring)
1106
1107	pLcm(f._poly, g._poly, m)
1108	p_Setm(m, self._ring)
1109	return new_MP(self,m)
1110
1111	def monomial_reduce(self, MPolynomial_libsingular f, G):
1112	"""
1113	Try to find a g in G where g.lm() divides f. If found (flt,g)
1114	is returned, (0,0) otherwise, where flt is f/g.lm().
1115
1116	It is assumed that G is iterable and contains ONLY elements in
1117	self.
1118
1119	INPUT:
1120	f -- monomial
1121	G -- list/set of mpolynomials
1122
1123	EXAMPLES:
1124	sage: P.<x,y,z> = QQ[]
1125	sage: f = x*y^2
1126	sage: G = [ 3/2x^3 + y^2 + 1/2, 1/4x*y + 2/7, 1/2 ]
1127	sage: P.monomial_reduce(f,G)
1128	(y, 1/4xy + 2/7)
1129
1130	TESTS:
1131	sage: P.<x,y,z> = QQ[]
1132	sage: f = x*y^2
1133	sage: G = [ 3/2x^3 + y^2 + 1/2, 1/4x*y + 2/7, 1/2 ]
1134
1135	sage: P.monomial_reduce(P(0),G)
1136	(0, 0)
1137
1138	sage: P.monomial_reduce(f,[P(0)])
1139	(0, 0)
1140
1141	"""
1142	cdef poly *m = f._poly
1143	cdef ring *r = self._ring
1144	cdef poly *flt
1145
1146	if not m:
1147	return f,f
1148
1149	for g in G:
1150	if PY_TYPE_CHECK(g, MPolynomial_libsingular) \
1151	and (<MPolynomial_libsingular>g) \
1152	and p_LmDivisibleBy((<MPolynomial_libsingular>g)._poly, m, r):
1153	flt = pDivide(f._poly, (<MPolynomial_libsingular>g)._poly)
1154	#p_SetCoeff(flt, n_Div( p_GetCoeff(f._poly, r) , p_GetCoeff((<MPolynomial_libsingular>g)._poly, r), r), r)
1155	p_SetCoeff(flt, n_Init(1, r), r)
1156	return new_MP(self,flt), g
1157	return self._zero,self._zero
1158
1159	def monomial_pairwise_prime(self, MPolynomial_libsingular g, MPolynomial_libsingular h):
1160	"""
1161	Return True if h and g are pairwise prime. Both are treated as monomials.
1162
1163	INPUT:
1164	h -- monomial
1165	g -- monomial
1166
1167	EXAMPLES:
1168	sage: P.<x,y,z> = QQ[]
1169	sage: P.monomial_pairwise_prime(x^2*z^3, y^4)
1170	True
1171
1172	sage: P.monomial_pairwise_prime(1/2x^3y^2, 3/4*y^3)
1173	False
1174
1175	TESTS:
1176	sage: Q.<x,y,z> = QQ[]
1177	sage: P.<x,y,z> = QQ[]
1178	sage: P.monomial_pairwise_prime(x^2*z^3, Q('y^4'))
1179	True
1180
1181	sage: P.monomial_pairwise_prime(1/2x^3y^2, Q(0))
1182	True
1183
1184	sage: P.monomial_pairwise_prime(P(1/2),x)
1185	False
1186
1187
1188	"""
1189	cdef int i
1190	cdef ring *r
1191	cdef poly p, q
1192
1193	if h._parent is not g._parent:
1194	g = (<MPolynomialRing_libsingular>h._parent)._coerce_c(g)
1195
1196	r = (<MPolynomialRing_libsingular>h._parent)._ring
1197	p = g._poly
1198	q = h._poly
1199
1200	if p == NULL:
1201	if q == NULL:
1202	return False #GCD(0,0) = 0
1203	else:
1204	return True #GCD(x,0) = 1
1205
1206	elif q == NULL:
1207	return True # GCD(0,x) = 1
1208
1209	elif p_IsConstant(p,r) or p_IsConstant(q,r): # assuming a base field
1210	return False
1211
1212	for i from 1 <= i <= r.N:
1213	if p_GetExp(p,i,r) and p_GetExp(q,i,r):
1214	return False
1215	return True
1216
1217
1218
1219	def monomial_all_divisors(self, MPolynomial_libsingular t):
1220	"""
1221	Return a list of all monomials that divide t, coefficients are
1222	ignored.
1223
1224	INPUT:
1225	t -- a monomial
1226
1227	OUTPUT:
1228	a list of monomials
1229
1230
1231	EXAMPLE:
1232	sage: P.<x,y,z> = QQ[]
1233	sage: P.monomial_all_divisors(x^2*z^3)
1234	[x, x^2, z, xz, x^2z, z^2, xz^2, x^2z^2, z^3, xz^3, x^2z^3]
1235
1236	ALGORITHM: addwithcarry idea by Toon Segers
1237	"""
1238
1239	M = list()
1240
1241	cdef ring *_ring = self._ring
1242	cdef poly *maxvector = t._poly
1243	cdef poly *tempvector = p_ISet(1, _ring)
1244
1245	pos = 1
1246
1247	while not p_ExpVectorEqual(tempvector, maxvector, _ring):
1248	tempvector = addwithcarry(tempvector, maxvector, pos, _ring)
1249	M.append(new_MP(self, p_Copy(tempvector,_ring)))
1250	return M
1251
1252
1253
1254	def unpickle_MPolynomialRing_libsingular(base_ring, names, term_order):
1255	"""
1256	inverse function for MPolynomialRing_libsingular.__reduce__
1257
1258	"""
1259	from sage.rings.polynomial.polynomial_ring_constructor import _get_from_cache
1260	key = (base_ring, tuple(names), len(names), False, term_order)
1261	R = _get_from_cache(key)
1262	if not R is None:
1263	return R
1264
1265	return MPolynomialRing_libsingular(base_ring, len(names), names, term_order)
1266
1267	cdef MPolynomial_libsingular new_MP(MPolynomialRing_libsingular parent, poly *juice):
1268	"""
1269	Construct a new MPolynomial_libsingular element
1270	"""
1271	cdef MPolynomial_libsingular p
1272	p = PY_NEW(MPolynomial_libsingular)
1273	p._parent = <ParentWithBase>parent
1274	p._poly = juice
1275	return p
1276
1277	cdef class MPolynomial_libsingular(sage.rings.polynomial.multi_polynomial.MPolynomial):
1278	"""
1279	A multivariate polynomial implemented using libSINGULAR.
1280	"""
1281	def __init__(self, MPolynomialRing_libsingular parent):
1282	"""
1283	Construct a zero element in parent.
1284	"""
1285	self._poly = NULL
1286	self._parent = <ParentWithBase>parent
1287
1288	def __dealloc__(self):
1289	# for some mysterious reason, various things may be NULL in some cases
1290	if self._parent is not <ParentWithBase>None and (<MPolynomialRing_libsingular>self._parent)._ring != NULL and self._poly != NULL:
1291	p_Delete(&self._poly, (<MPolynomialRing_libsingular>self._parent)._ring)
1292
1293	def __call__(self, x, *kwds):
1294	r"""
1295	Evaluate this multi-variate polynomial at $x$, where $x$ is
1296	either the tuple of values to substitute in, or one can use
1297	functional notation $f(a_0,a_1,a_2, \ldots)$ to evaluate $f$
1298	with the ith variable replaced by $a_i$.
1299
1300	INPUT:
1301	x -- a list of elements in self.parent()
1302	or **kwds -- a dictionary of variable-name:value pairs.
1303
1304	EXAMPLES:
1305	sage: P.<x,y,z> = QQ[]
1306	sage: f = 3/2x^2y + 1/7 * y^2 + 13/27
1307	sage: f(0,0,0)
1308	13/27
1309
1310	sage: f(1,1,1)
1311	803/378
1312	sage: 3/2 + 1/7 + 13/27
1313	803/378
1314
1315	sage: f(45/2,19/3,1)
1316	7281167/1512
1317
1318	TESTS:
1319	sage: P.<x,y,z> = QQ[]
1320	sage: P(0)(1,2,3)
1321	0
1322	sage: P(3/2)(1,2,3)
1323	3/2
1324
1325	sage: R.<a,b,y> = QQ[]
1326	sage: f = ay^3 + by - 3ab*y
1327	sage: f(a=5,b=3,y=10)
1328	4580
1329	sage: f(5,3,10)
1330	4580
1331	"""
1332	if len(kwds) > 0:
1333	f = self.subs(**kwds)
1334	if len(x) > 0:
1335	return f(*x)
1336	else:
1337	return f
1338
1339	cdef int l = len(x)
1340	cdef MPolynomialRing_libsingular parent = (<MPolynomialRing_libsingular>self._parent)
1341	cdef ring *_ring = parent._ring
1342
1343	cdef poly *_p
1344
1345	if l == 1 and (PY_TYPE_CHECK(x, tuple) or PY_TYPE_CHECK(x, list)):
1346	x = x[0]
1347	l = len(x)
1348
1349	if l != parent._ring.N:
1350	raise TypeError, "number of arguments does not match number of variables in parent"
1351
1352	try:
1353	x = [parent._coerce_c(e) for e in x]
1354	except TypeError:
1355	# give up, evaluate functional
1356	y = parent.base_ring()(0)
1357	for (m,c) in self.dict().iteritems():
1358	y += cmul([ x[i]*m[i] for i in m.nonzero_positions()])
1359	return y
1360
1361	cdef ideal *to_id = idInit(l,1)
1362	for i from 0 <= i < l:
1363	to_id.m[i]= p_Copy( (<MPolynomial_libsingular>x[i])._poly, _ring)
1364
1365	cdef ideal *from_id=idInit(1,1)
1366	from_id.m[0] = self._poly
1367
1368	cdef ideal *res_id = fast_map(from_id, _ring, to_id, _ring)
1369	cdef poly *res = res_id.m[0]
1370
1371	from_id.m[0] = NULL
1372	res_id.m[0] = NULL
1373
1374	id_Delete(&to_id, _ring)
1375	id_Delete(&from_id, _ring)
1376	id_Delete(&res_id, _ring)
1377	return new_MP(parent, res)
1378
1379	def __richcmp__(left, right, int op):
1380	return (<Element>left)._richcmp(right, op)
1381
1382	cdef int _cmp_c_impl(left, Element right) except -2:
1383	"""
1384	Compare left and right and return -1, 0, and 1 for <,==, and > respectively.
1385
1386	EXAMPLES:
1387	sage: P.<x,y,z> = MPolynomialRing(QQ,3, order='degrevlex')
1388	sage: x == x
1389	True
1390
1391	sage: x > y
1392	True
1393	sage: y^2 > x
1394	True
1395
1396	sage: (2/3x^2 + 1/2y + 3) > (2/3x^2 + 1/4y + 10)
1397	True
1398
1399	TESTS:
1400	sage: P.<x,y,z> = MPolynomialRing(QQ,3, order='degrevlex')
1401	sage: x > P(0)
1402	True
1403
1404	sage: P(0) == P(0)
1405	True
1406
1407	sage: P(0) < P(1)
1408	True
1409
1410	sage: x > P(1)
1411	True
1412
1413	sage: 1/2x < 3/4x
1414	True
1415
1416	sage: (x+1) > x
1417	True
1418
1419	sage: f = 3/4x^2y + 1/2*x + 2/7
1420	sage: f > f
1421	False
1422	sage: f < f
1423	False
1424	sage: f == f
1425	True
1426
1427	sage: P.<x,y,z> = MPolynomialRing(GF(127),3, order='degrevlex')
1428	sage: (66x^2 + 23) > (66x^2 + 2)
1429	True
1430
1431
1432	"""
1433	cdef ring *r
1434	cdef poly p, q
1435	cdef number *h
1436	cdef int ret = 0
1437
1438	r = (<MPolynomialRing_libsingular>left._parent)._ring
1439	if(r != currRing): rChangeCurrRing(r)
1440	p = (<MPolynomial_libsingular>left)._poly
1441	q = (<MPolynomial_libsingular>right)._poly
1442
1443	# handle special cases first (slight slowdown, as special
1444	# cases are - well - special
1445	if p==NULL:
1446	if q==NULL:
1447	# compare 0, 0
1448	return 0
1449	elif p_IsConstant(q,r):
1450	# compare 0, const
1451	return 1-2*n_GreaterZero(p_GetCoeff(q,r), r) # -1: <, 1: > #
1452	elif q==NULL:
1453	if p_IsConstant(p,r):
1454	# compare const, 0
1455	return -1+2*n_GreaterZero(p_GetCoeff(p,r), r) # -1: <, 1: >
1456	#else
1457
1458	while ret==0 and p!=NULL and q!=NULL:
1459	ret = p_Cmp( p, q, r)
1460
1461	if ret==0:
1462	h = n_Sub(p_GetCoeff(p, r),p_GetCoeff(q, r), r)
1463	# compare coeffs
1464	ret = -1+n_IsZero(h, r)+2*n_GreaterZero(h, r) # -1: <, 0:==, 1: >
1465	n_Delete(&h, r)
1466	p = pNext(p)
1467	q = pNext(q)
1468
1469	if ret==0:
1470	if p==NULL and q != NULL:
1471	ret = -1
1472	elif p!=NULL and q==NULL:
1473	ret = 1
1474
1475	return ret
1476
1477	cdef ModuleElement _add_c_impl( left, ModuleElement right):
1478	"""
1479	Add left and right.
1480
1481	EXAMPLE:
1482	sage: P.<x,y,z>=MPolynomialRing(QQ,3)
1483	sage: 3/2x + 1/2y + 1
1484	3/2x + 1/2y + 1
1485
1486	"""
1487	cdef MPolynomial_libsingular res
1488
1489	cdef poly _l, _r, *_p
1490	cdef ring *_ring
1491
1492	_ring = (<MPolynomialRing_libsingular>left._parent)._ring
1493
1494	_l = p_Copy(left._poly, _ring)
1495	_r = p_Copy((<MPolynomial_libsingular>right)._poly, _ring)
1496
1497	if(_ring != currRing): rChangeCurrRing(_ring)
1498	_p= p_Add_q(_l, _r, _ring)
1499
1500	return new_MP((<MPolynomialRing_libsingular>left._parent),_p)
1501
1502	cdef ModuleElement _sub_c_impl( left, ModuleElement right):
1503	"""
1504	Subtract left and right.
1505
1506	EXAMPLE:
1507	sage: P.<x,y,z>=MPolynomialRing(QQ,3)
1508	sage: 3/2x - 1/2y - 1
1509	3/2x - 1/2y - 1
1510
1511	"""
1512	cdef MPolynomial_libsingular res
1513
1514	cdef poly _l, _r, *_p
1515	cdef ring *_ring
1516
1517	_ring = (<MPolynomialRing_libsingular>left._parent)._ring
1518
1519	_l = p_Copy(left._poly, _ring)
1520	_r = p_Copy((<MPolynomial_libsingular>right)._poly, _ring)
1521
1522	if(_ring != currRing): rChangeCurrRing(_ring)
1523	_p= p_Add_q(_l, p_Neg(_r, _ring), _ring)
1524
1525	return new_MP((<MPolynomialRing_libsingular>left._parent),_p)
1526
1527
1528	cdef ModuleElement _rmul_c_impl(self, RingElement left):
1529	"""
1530	Multiply self with a base ring element.
1531
1532	EXAMPLE:
1533	sage: P.<x,y,z>=MPolynomialRing(QQ,3)
1534	sage: 3/2*x
1535	3/2*x
1536	"""
1537
1538	cdef number *_n
1539	cdef ring *_ring
1540	cdef poly *_p
1541
1542	_ring = (<MPolynomialRing_libsingular>self._parent)._ring
1543
1544	if(_ring != currRing): rChangeCurrRing(_ring)
1545
1546	if not left:
1547	return (<MPolynomialRing_libsingular>self._parent)._zero
1548
1549	_n = co.sa2si(left,_ring)
1550
1551	_p = pp_Mult_nn(self._poly,_n,_ring)
1552	n_Delete(&_n, _ring)
1553	return new_MP((<MPolynomialRing_libsingular>self._parent),_p)
1554
1555	cdef RingElement _mul_c_impl(left, RingElement right):
1556	"""
1557	Multiply left and right.
1558
1559	EXAMPLE:
1560	sage: P.<x,y,z>=MPolynomialRing(QQ,3)
1561	sage: (3/2x - 1/2y - 1) * (3/2x + 1/2y + 1)
1562	9/4x^2 - 1/4y^2 - y - 1
1563	"""
1564	cdef poly _l, _r, *_p
1565	cdef ring *_ring
1566
1567	_ring = (<MPolynomialRing_libsingular>left._parent)._ring
1568
1569	if(_ring != currRing): rChangeCurrRing(_ring)
1570	_p = pp_Mult_qq(left._poly, (<MPolynomial_libsingular>right)._poly, _ring)
1571	return new_MP(left._parent,_p)
1572
1573	cdef RingElement _div_c_impl(left, RingElement right):
1574	"""
1575	Divide left by right
1576
1577	EXAMPLES:
1578	sage: R.<x,y>=MPolynomialRing(QQ,2)
1579	sage: f = (x + y)/3
1580	sage: f.parent()
1581	Polynomial Ring in x, y over Rational Field
1582
1583	Note that / is still a constructor for elements of the
1584	fraction field in all cases as long as both arguments have the
1585	same parent.
1586
1587	sage: R.<x,y>=PolynomialRing(QQ,2)
1588	sage: R.<x,y>=MPolynomialRing(QQ,2)
1589	sage: f = x^3 + y
1590	sage: g = x
1591	sage: h = f/g; h
1592	(x^3 + y)/x
1593	sage: h.parent()
1594	Fraction Field of Polynomial Ring in x, y over Rational Field
1595
1596	TESTS:
1597	sage: R.<x,y>=MPolynomialRing(QQ,2)
1598	sage: x/0
1599	Traceback (most recent call last):
1600	...
1601	ZeroDivisionError
1602
1603	"""
1604	cdef poly *p
1605	cdef ring *r
1606	cdef number *n
1607	if (<MPolynomial_libsingular>right).is_constant_c():
1608
1609	p = (<MPolynomial_libsingular>right)._poly
1610	if p == NULL:
1611	raise ZeroDivisionError
1612	r = (<MPolynomialRing_libsingular>(<MPolynomial_libsingular>left)._parent)._ring
1613	n = p_GetCoeff(p, r)
1614	n = nInvers(n)
1615	p = pp_Mult_nn(left._poly, n, r)
1616	n_Delete(&n,r)
1617	return new_MP(left._parent, p)
1618	else:
1619	return (<MPolynomialRing_libsingular>left._parent).fraction_field()(left,right)
1620
1621	def __pow__(MPolynomial_libsingular self,int exp,ignored):
1622	"""
1623	Return self^(exp).
1624
1625	EXAMPLE:
1626	sage: R.<x,y>=MPolynomialRing(QQ,2)
1627	sage: f = x^3 + y
1628	sage: f^2
1629	x^6 + 2x^3y + y^2
1630	sage: g = f^(-1); g
1631	1/(x^3 + y)
1632	sage: type(g)
1633	<class 'sage.rings.fraction_field_element.FractionFieldElement'>
1634	"""
1635	cdef ring *_ring
1636	_ring = (<MPolynomialRing_libsingular>self._parent)._ring
1637
1638	cdef poly *_p
1639	cdef int _exp
1640
1641	_exp = exp
1642
1643	if _exp < 0:
1644	return 1/(self**(-_exp))
1645	if _exp > 65535:
1646	raise TypeError, "exponent is too large, max. is 65535"
1647
1648	if(_ring != currRing): rChangeCurrRing(_ring)
1649
1650	_p = pPower( p_Copy(self._poly,_ring),_exp)
1651
1652	return new_MP((<MPolynomialRing_libsingular>self._parent),_p)
1653
1654
1655	def __neg__(self):
1656	"""
1657	Return -self.
1658
1659	EXAMPLE:
1660	sage: R.<x,y>=MPolynomialRing(QQ,2)
1661	sage: f = x^3 + y
1662	sage: -f
1663	-x^3 - y
1664	"""
1665	cdef ring *_ring
1666	_ring = (<MPolynomialRing_libsingular>self._parent)._ring
1667	if(_ring != currRing): rChangeCurrRing(_ring)
1668
1669	return new_MP((<MPolynomialRing_libsingular>self._parent),\
1670	p_Neg(p_Copy(self._poly,_ring),_ring))
1671
1672	def _repr_(self):
1673	s = self._repr_short_c()
1674	s = s.replace("+"," + ").replace("-"," - ")
1675	if s.startswith(" - "):
1676	return "-" + s[3:]
1677	else:
1678	return s
1679
1680	def _repr_short(self):
1681	"""
1682	This is a faster but less pretty way to print polynomials. If available
1683	it uses the short SINGULAR notation.
1684	"""
1685	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
1686	if _ring.CanShortOut:
1687	_ring.ShortOut = 1
1688	s = self._repr_short_c()
1689	_ring.ShortOut = 0
1690	else:
1691	s = self._repr_short_c()
1692	return s
1693
1694	cdef _repr_short_c(self):
1695	"""
1696	Raw SINGULAR printing.
1697	"""
1698	rChangeCurrRing((<MPolynomialRing_libsingular>self._parent)._ring)
1699	s = p_String(self._poly, (<MPolynomialRing_libsingular>self._parent)._ring, (<MPolynomialRing_libsingular>self._parent)._ring)
1700	return s
1701
1702	def _latex_(self):
1703	r"""
1704	Return a polynomial latex representation of self.
1705
1706	EXAMPLE:
1707	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
1708	sage: f = - 1x^2y - 25/27 * y^3 - z^2
1709	sage: latex(f)
1710	- x^{2}y - \frac{25}{27} y^{3} - z^{2}
1711
1712	"""
1713	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
1714	cdef int n = _ring.N
1715	cdef int j, e
1716	cdef poly *p = self._poly
1717	poly = ""
1718	gens = self.parent().latex_variable_names()
1719	base = self.parent().base()
1720
1721	while p:
1722	sign_switch = False
1723
1724	# First determine the multinomial:
1725	multi = ""
1726	for j from 1 <= j <= n:
1727	e = p_GetExp(p, j, _ring)
1728	if e > 0:
1729	multi += gens[j-1]
1730	if e > 1:
1731	multi += "^{%d}"%e
1732
1733	# Next determine coefficient of multinomial
1734	c = co.si2sa( p_GetCoeff(p, _ring), _ring, base)
1735	if len(multi) == 0:
1736	multi = latex(c)
1737	elif c != 1:
1738	if c == -1:
1739	if len(poly) > 0:
1740	sign_switch = True
1741	else:
1742	multi = "- %s"%(multi)
1743	else:
1744	multi = "%s %s"%(latex(c),multi)
1745
1746	# Now add on coefficiented multinomials
1747	if len(poly) > 0:
1748	if sign_switch:
1749	poly = poly + " - "
1750	else:
1751	poly = poly + " + "
1752	poly = poly + multi
1753
1754	p = pNext(p)
1755
1756	poly = poly.lstrip().rstrip()
1757	poly = poly.replace("+ -","- ")
1758
1759	if len(poly) == 0:
1760	return "0"
1761	return poly
1762
1763	def _macaulay2_(self, macaulay2=macaulay2):
1764	"""
1765	Return corresponding Macaulay2 polynomial.
1766
1767	WARNING: Two identical rings are not canonically isomorphic in
1768	M2, so we require the user to explicitly set the ring, since
1769	there is no way to know if the ring has been set or not, and
1770	setting it twice screws everything up.
1771
1772	EXAMPLES:
1773	sage: R.<x,y> = PolynomialRing(GF(7), 2) # optional
1774	sage: f = (x^3 + 2y^2x)^7; f # optional
1775	x^21 + 2x^7y^14
1776
1777	Always call the Macaulay2 ring conversion on the parent polynomial
1778	ring before converting a copy of elements to Macaulay2:
1779	sage: macaulay2(R) # optional
1780	ZZ/7 [x, y, MonomialOrder => GRevLex, MonomialSize => 16]
1781	sage: h = f._macaulay2_(); h # optional
1782	x^21+2x^7y^14
1783	sage: k = (x+y)._macaulay2_() # optional
1784	sage: k + h # optional
1785	x^21+2x^7y^14+x+y
1786	sage: R(h) # optional
1787	x^21 + 2x^7y^14
1788	sage: R(h^20) == f^20 # optional
1789	True
1790	"""
1791	try:
1792	if self.__macaulay2[macaulay2].parent() is macaulay2:
1793	return self.__macaulay2[macaulay2]
1794	except (TypeError, AttributeError):
1795	self.__macaulay2 = {}
1796	except KeyError:
1797	pass
1798	#self.parent()._macaulay2_set_ring(macaulay2)
1799	z = macaulay2(repr(self))
1800	self.__macaulay2[macaulay2] = z
1801	return z
1802
1803	def _repr_with_changed_varnames(self, varnames):
1804	"""
1805	Return string representing self but change the variable names
1806	to varnames.
1807
1808	EXAMPLE:
1809	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
1810	sage: f = - 1x^2y - 25/27 * y^3 - z^2
1811	sage: print f._repr_with_changed_varnames(['FOO', 'BAR', 'FOOBAR'])
1812	-FOO^2BAR - 25/27BAR^3 - FOOBAR^2
1813
1814	"""
1815	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
1816	cdef char **_names
1817	cdef char **_orig_names
1818	cdef char *_name
1819	cdef int i
1820
1821	if len(varnames) != _ring.N:
1822	raise TypeError, "len(varnames) doesn't equal self.parent().ngens()"
1823
1824
1825	_names = <char*>sage_malloc(sizeof(char)*_ring.N)
1826	for i from 0 <= i < _ring.N:
1827	_name = varnames[i]
1828	_names[i] = strdup(_name)
1829
1830	_orig_names = _ring.names
1831	_ring.names = _names
1832	s = str(self)
1833	_ring.names = _orig_names
1834
1835	for i from 0 <= i < _ring.N:
1836	free(_names[i]) # strdup() --> free()
1837	sage_free(_names)
1838
1839	return s
1840
1841	def degree(self, MPolynomial_libsingular x=None):
1842	"""
1843	Return the maximal degree of self in x, where x must be one of the
1844	generators for the parent of self.
1845
1846	INPUT:
1847	x -- multivariate polynmial (a generator of the parent of self)
1848	If x is not specified (or is None), return the total degree,
1849	which is the maximum degree of any monomial.
1850
1851	OUTPUT:
1852	integer
1853
1854	EXAMPLE:
1855	sage: R.<x, y> = MPolynomialRing(QQ, 2)
1856	sage: f = y^2 - x^9 - x
1857	sage: f.degree(x)
1858	9
1859	sage: f.degree(y)
1860	2
1861	sage: (y^10x - 7x^2y^5 + 5x^3).degree(x)
1862	3
1863	sage: (y^10x - 7x^2y^5 + 5x^3).degree(y)
1864	10
1865
1866	TESTS:
1867	sage: P.<x, y> = MPolynomialRing(QQ, 2)
1868	sage: P(0).degree(x)
1869	0
1870	sage: P(1).degree(x)
1871	0
1872
1873	"""
1874	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
1875	cdef poly *p = self._poly
1876	cdef int deg, _deg
1877
1878	deg = 0
1879
1880	if not x:
1881	return self.total_degree()
1882
1883	# TODO: we can do this faster
1884	if not x in self._parent.gens():
1885	raise TypeError, "x must be one of the generators of the parent."
1886	for i from 1 <= i <= r.N:
1887	if p_GetExp(x._poly, i, r):
1888	break
1889	while p:
1890	_deg = p_GetExp(p,i,r)
1891	if _deg > deg:
1892	deg = _deg
1893	p = pNext(p)
1894
1895	return deg
1896
1897	def newton_polytope(self):
1898	"""
1899	Return the Newton polytope of this polynomial.
1900
1901	You should have the optional polymake package installed.
1902
1903	EXAMPLES:
1904	sage: R.<x,y> = MPolynomialRing(QQ,2)
1905	sage: f = 1 + x*y + x^3 + y^3
1906	sage: P = f.newton_polytope()
1907	sage: P
1908	Convex hull of points [[1, 0, 0], [1, 0, 3], [1, 1, 1], [1, 3, 0]]
1909	sage: P.facets()
1910	[(0, 1, 0), (3, -1, -1), (0, 0, 1)]
1911	sage: P.is_simple()
1912	True
1913
1914	TESTS:
1915	sage: R.<x,y> = MPolynomialRing(QQ,2)
1916	sage: R(0).newton_polytope()
1917	Convex hull of points []
1918	sage: R(1).newton_polytope()
1919	Convex hull of points [[1, 0, 0]]
1920
1921	"""
1922	from sage.geometry.all import polymake
1923	e = self.exponents()
1924	a = [[1] + list(v) for v in e]
1925	P = polymake.convex_hull(a)
1926	return P
1927
1928	def total_degree(self):
1929	"""
1930	Return the total degree of self, which is the maximum degree
1931	of all monomials in self.
1932
1933	EXAMPLES:
1934	sage: R.<x,y,z> = MPolynomialRing(QQ, 3)
1935	sage: f=2xy^3*z^2
1936	sage: f.total_degree()
1937	6
1938	sage: f=4x^2y^2*z^3
1939	sage: f.total_degree()
1940	7
1941	sage: f=99x^6y^3*z^9
1942	sage: f.total_degree()
1943	18
1944	sage: f=xy^3z^6+3*x^2
1945	sage: f.total_degree()
1946	10
1947	sage: f=z^3+8x^4y^5*z
1948	sage: f.total_degree()
1949	10
1950	sage: f=z^9+10x^4+y^8x^2
1951	sage: f.total_degree()
1952	10
1953
1954	TESTS:
1955	sage: R.<x,y,z> = MPolynomialRing(QQ, 3)
1956	sage: R(0).total_degree()
1957	0
1958	sage: R(1).total_degree()
1959	0
1960	"""
1961	cdef poly *p = self._poly
1962	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
1963	cdef int l
1964	if self._poly == NULL:
1965	return 0
1966	if(r != currRing): rChangeCurrRing(r)
1967	return pLDeg(p,&l,r)
1968
1969	def monomial_coefficient(self, MPolynomial_libsingular mon):
1970	"""
1971	Return the coefficient of the monomial mon in self, where mon
1972	must have the same parent as self.
1973
1974	INPUT:
1975	mon -- a monomial
1976
1977	OUTPUT:
1978	ring element
1979
1980	EXAMPLE:
1981	sage: P.<x,y> = MPolynomialRing(QQ, 2)
1982
1983	The coefficient returned is an element of the base ring of self; in
1984	this case, QQ.
1985	sage: f = 2 * x * y
1986	sage: c = f.monomial_coefficient(x*y); c
1987	2
1988	sage: c in QQ
1989	True
1990
1991	sage: f = y^2 - x^9 - 7x + 5x*y
1992	sage: f.monomial_coefficient(y^2)
1993	1
1994	sage: f.monomial_coefficient(x*y)
1995	5
1996	sage: f.monomial_coefficient(x^9)
1997	-1
1998	sage: f.monomial_coefficient(x^10)
1999	0
2000	"""
2001	cdef poly *p = self._poly
2002	cdef poly *m = mon._poly
2003	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2004
2005	if not mon._parent is self._parent:
2006	raise TypeError, "mon must have same parent as self"
2007
2008	while(p):
2009	if p_ExpVectorEqual(p, m, r) == 1:
2010	return co.si2sa(p_GetCoeff(p, r), r, (<MPolynomialRing_libsingular>self._parent)._base)
2011	p = pNext(p)
2012
2013	return (<MPolynomialRing_libsingular>self._parent)._base._zero_element
2014
2015	def dict(self):
2016	"""
2017	Return a dictionary representing self. This dictionary is in
2018	the same format as the generic MPolynomial: The dictionary
2019	consists of ETuple:coefficient pairs.
2020
2021	EXAMPLE:
2022	sage: R.<x,y,z> = MPolynomialRing(QQ, 3)
2023	sage: f=2xy^3z^2 + 1/7x^2 + 2/3
2024	sage: f.dict()
2025	{(2, 0, 0): 1/7, (0, 0, 0): 2/3, (1, 3, 2): 2}
2026	"""
2027	cdef poly *p
2028	cdef ring *r
2029	cdef int n
2030	cdef int v
2031	r = (<MPolynomialRing_libsingular>self._parent)._ring
2032	if r!=currRing: rChangeCurrRing(r)
2033	base = (<MPolynomialRing_libsingular>self._parent)._base
2034	p = self._poly
2035	pd = dict()
2036	while p:
2037	d = dict()
2038	for v from 1 <= v <= r.N:
2039	n = p_GetExp(p,v,r)
2040	if n!=0:
2041	d[v-1] = n
2042
2043	pd[ETuple(d,r.N)] = co.si2sa(p_GetCoeff(p, r), r, base)
2044
2045	p = pNext(p)
2046	return pd
2047
2048	def __getitem__(self,x):
2049	"""
2050	same as self.monomial_coefficent but for exponent vectors.
2051
2052	INPUT:
2053	x -- a tuple or, in case of a single-variable MPolynomial
2054	ring x can also be an integer.
2055
2056	EXAMPLES:
2057	sage: R.<x, y> = MPolynomialRing(QQ, 2)
2058	sage: f = -10x^3y + 17xy
2059	sage: f[3,1]
2060	-10
2061	sage: f[1,1]
2062	17
2063	sage: f[0,1]
2064	0
2065
2066	sage: R.<x> = MPolynomialRing(GF(7),1); R
2067	Polynomial Ring in x over Finite Field of size 7
2068	sage: f = 5*x^2 + 3; f
2069	-2*x^2 + 3
2070	sage: f[2]
2071	5
2072	"""
2073
2074	cdef poly *m
2075	cdef poly *p = self._poly
2076	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2077	cdef int i
2078
2079	if PY_TYPE_CHECK(x, MPolynomial_libsingular):
2080	return self.monomial_coefficient(x)
2081	if not PY_TYPE_CHECK(x, tuple):
2082	try:
2083	x = tuple(x)
2084	except TypeError:
2085	x = (x,)
2086
2087	if len(x) != (<MPolynomialRing_libsingular>self._parent).__ngens:
2088	raise TypeError, "x must have length self.ngens()"
2089
2090	m = p_ISet(1,r)
2091	i = 1
2092	for e in x:
2093	p_SetExp(m, i, int(e), r)
2094	i += 1
2095	p_Setm(m, r)
2096
2097	while(p):
2098	if p_ExpVectorEqual(p, m, r) == 1:
2099	p_Delete(&m,r)
2100	return co.si2sa(p_GetCoeff(p, r), r, (<MPolynomialRing_libsingular>self._parent)._base)
2101	p = pNext(p)
2102
2103	p_Delete(&m,r)
2104	return (<MPolynomialRing_libsingular>self._parent)._base._zero_element
2105
2106	def coefficient(self, MPolynomial_libsingular mon):
2107	"""
2108	Return the coefficient of mon in self, where mon must have the
2109	same parent as self. The coefficient is defined as follows.
2110	If f is this polynomial, then the coefficient is the sum T/mon
2111	where the sum is over terms T in f that are exactly divisible
2112	by mon.
2113
2114	A monomial m(x,y) 'exactly divides' f(x,y) if m(x,y)\|f(x,y)
2115	and neither xm(x,y) nor ym(x,y) divides f(x,y).
2116
2117	INPUT:
2118	mon -- a monomial
2119
2120	OUTPUT:
2121	element of the parent of self
2122
2123	EXAMPLES:
2124	sage: P.<x,y> = MPolynomialRing(QQ, 2)
2125
2126	The coefficient returned is an element of the parent of self; in
2127	this case, QQ[x, y].
2128
2129	sage: f = 2 * x * y
2130	sage: c = f.coefficient(x*y); c
2131	2
2132	sage: c.parent()
2133	Polynomial Ring in x, y over Rational Field
2134	sage: c in P
2135	True
2136
2137	sage: f = y^2 - x^9 - 7x + 5x*y
2138	sage: f.coefficient(y)
2139	5*x
2140	sage: f = y - x^9y - 7x + 5xy
2141	sage: f.coefficient(y)
2142	-x^9 + 5*x + 1
2143
2144	sage: f = y^2 - x^9 - 7xy^2 + 5xy
2145	sage: f.coefficient(x*y)
2146	5
2147	sage: f.coefficient(x*y^2)
2148	-7
2149	sage: f.coefficient(y)
2150	5*x
2151	sage: f.coefficient(x)
2152	-7y^2 + 5y
2153
2154	The coefficient of 1 is also an element of the multivariate
2155	polynomial ring:
2156
2157	sage: R.<x,y> = MPolynomialRing(GF(389),2)
2158	sage: parent(R(x*y+5).coefficient(R(1)))
2159	Polynomial Ring in x, y over Finite Field of size 389
2160	"""
2161	cdef poly *p = self._poly
2162	cdef poly *m = mon._poly
2163	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2164	cdef poly *res = p_ISet(0,r)
2165	cdef poly *t
2166	cdef int exactly_divisible, i
2167	if(r != currRing): rChangeCurrRing(r)
2168
2169	if not mon._parent is self._parent:
2170	raise TypeError, "mon must have same parent as self"
2171
2172	while(p):
2173	exactly_divisible = 1
2174	for i from 1 <= i <= r.N:
2175	if (p_GetExp(m,i,r) != 0) and (p_GetExp(p,i,r) != p_GetExp(m,i,r)):
2176	exactly_divisible = 0
2177	break
2178	if exactly_divisible:
2179	t = pDivide(p,m)
2180	p_SetCoeff(t, n_Div( p_GetCoeff(p, r) , p_GetCoeff(m, r), r), r)
2181	res = p_Add_q(res, t , r )
2182	p = pNext(p)
2183
2184	return new_MP(self._parent, res)
2185
2186	def exponents(self):
2187	"""
2188	Return the exponents of the monomials appearing in self.
2189
2190	EXAMPLES:
2191	sage: R.<a,b,c> = PolynomialRing(QQ, 3)
2192	sage: R.<a,b,c> = MPolynomialRing(QQ, 3)
2193	sage: f = a^3 + b + 2*b^2
2194	sage: f.exponents()
2195	[(3, 0, 0), (0, 2, 0), (0, 1, 0)]
2196	"""
2197	cdef poly *p
2198	cdef ring *r
2199	cdef int v
2200	r = (<MPolynomialRing_libsingular>self._parent)._ring
2201
2202	p = self._poly
2203
2204	pl = list()
2205	while p:
2206	ml = list()
2207	for v from 1 <= v <= r.N:
2208	ml.append(p_GetExp(p,v,r))
2209	pl.append(ETuple(ml))
2210
2211	p = pNext(p)
2212	return pl
2213
2214	def is_unit(self):
2215	"""
2216	Return True if self is a unit.
2217
2218	EXAMPLES:
2219	sage: R.<x,y> = PolynomialRing(QQ, 2)
2220	sage: R.<x,y> = MPolynomialRing(QQ, 2)
2221	sage: (x+y).is_unit()
2222	False
2223	sage: R(0).is_unit()
2224	False
2225	sage: R(-1).is_unit()
2226	True
2227	sage: R(-1 + x).is_unit()
2228	False
2229	sage: R(2).is_unit()
2230	True
2231	"""
2232	return bool(p_IsUnit(self._poly, (<MPolynomialRing_libsingular>self._parent)._ring))
2233
2234	def inverse_of_unit(self):
2235	"""
2236	Return the inverse of self if self is a unit.
2237
2238	EXAMPLES:
2239
2240	sage: R.<x,y> = MPolynomialRing(QQ, 2)
2241	"""
2242	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
2243	if(_ring != currRing): rChangeCurrRing(_ring)
2244
2245	if not p_IsUnit(self._poly, _ring):
2246	raise ArithmeticError, "is not a unit"
2247	else:
2248	return new_MP(self._parent,pInvers(0,self._poly,NULL))
2249
2250	def is_homogeneous(self):
2251	"""
2252	Return True if self is a homogeneous polynomial.
2253
2254	EXAMPLES:
2255	sage: P.<x,y> = PolynomialRing(RationalField(), 2)
2256	sage: (x+y).is_homogeneous()
2257	True
2258	sage: (x.parent()(0)).is_homogeneous()
2259	True
2260	sage: (x+y^2).is_homogeneous()
2261	False
2262	sage: (x^2 + y^2).is_homogeneous()
2263	True
2264	sage: (x^2 + y^2*x).is_homogeneous()
2265	False
2266	sage: (x^2y + y^2x).is_homogeneous()
2267	True
2268
2269	"""
2270	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
2271	if(_ring != currRing): rChangeCurrRing(_ring)
2272	return bool(pIsHomogeneous(self._poly))
2273
2274	def homogenize(self, var='h'):
2275	"""
2276	Return self is self is homogeneous. Otherwise return a
2277	homogeneous polynomial. If a string is given, return a
2278	polynomial in one more variable such that setting that
2279	variable equal to 1 yields self. This variable is added to the
2280	end of the variables. If either a variable in self.parent() or
2281	an index is given, this variable is used to homogenize the
2282	polynomial.
2283
2284	INPUT:
2285	var -- either a string (default: 'h'); a variable name for the new variable
2286	to be added in when homogenizing or a variable/index to specify the existing
2287	variable to be used.
2288
2289	OUTPUT:
2290	a multivariate polynomial
2291
2292	EXAMPLES:
2293	sage: P.<x,y> = PolynomialRing(QQ,2)
2294	sage: P.<x,y> = MPolynomialRing(QQ,2)
2295	sage: f = x^2 + y + 1 + 5xy^10
2296	sage: g = f.homogenize('z'); g
2297	5xy^10 + x^2z^9 + yz^10 + z^11
2298	sage: g.parent()
2299	Polynomial Ring in x, y, z over Rational Field
2300	sage: f.homogenize(x)
2301	2x^11 + x^10y + 5xy^10
2302
2303	"""
2304	cdef MPolynomialRing_libsingular parent = <MPolynomialRing_libsingular>self._parent
2305	cdef MPolynomial_libsingular f
2306
2307	if self.is_homogeneous():
2308	return self
2309
2310	if PY_TYPE_CHECK(var, MPolynomial_libsingular):
2311	if (<MPolynomial_libsingular>var)._parent is self._parent:
2312	var = var._variable_indices_()
2313	if len(var) == 1:
2314	var = var[0]
2315	else:
2316	raise TypeError, "parameter var must be single variable"
2317
2318	if PY_TYPE_CHECK(var,str):
2319	names = [str(e) for e in parent.gens()] + [var]
2320	P = MPolynomialRing_libsingular(parent.base(),parent.ngens()+1, names, order=parent.term_order())
2321	f = P(str(self))
2322	return new_MP(P, pHomogen(f._poly,len(names)))
2323	elif PY_TYPE_CHECK(var,int) or PY_TYPE_CHECK(var,Integer):
2324	if var < parent._ring.N:
2325	return new_MP(parent, pHomogen(p_Copy(self._poly, parent._ring), var+1))
2326	else:
2327	raise TypeError, "var must be < self.parent().ngens()"
2328	else:
2329	raise TypeError, "parameter var not understood"
2330
2331	def is_monomial(self):
2332	return not self._poly.next
2333
2334	def subs(self, fixed=None, **kw):
2335	"""
2336	Fixes some given variables in a given multivariate polynomial and
2337	returns the changed multivariate polynomials. The polynomial
2338	itself is not affected. The variable,value pairs for fixing are
2339	to be provided as dictionary of the form {variable:value}.
2340
2341	This is a special case of evaluating the polynomial with some of
2342	the variables constants and the others the original variables, but
2343	should be much faster if only few variables are to be fixed.
2344
2345	INPUT:
2346	fixed -- (optional) dict with variable:value pairs
2347	**kw -- names parameters
2348
2349	OUTPUT:
2350	new MPolynomial
2351
2352	EXAMPLES:
2353	sage: R.<x,y> = QQ[]
2354	sage: f = x^2 + y + x^2*y^2 + 5
2355	sage: f(5,y)
2356	25*y^2 + y + 30
2357	sage: f.subs({x:5})
2358	25*y^2 + y + 30
2359	sage: f.subs(x=5)
2360	25*y^2 + y + 30
2361
2362	TESTS:
2363	sage: P.<x,y,z> = QQ[]
2364	sage: f = y
2365	sage: f.subs({y:x}).subs({x:z})
2366	z
2367
2368	NOTE: The evaluation is performed by evalutating every
2369	variable:value pair separately. This has side effects if
2370	e.g. x=y, y=z is provided. If x=y is evaluated first, all x
2371	variables will be replaced by z eventually.
2372
2373	"""
2374	cdef int mi, i
2375
2376	cdef MPolynomialRing_libsingular parent = <MPolynomialRing_libsingular>self._parent
2377	cdef ring *_ring = parent._ring
2378
2379	if(_ring != currRing): rChangeCurrRing(_ring)
2380
2381	cdef poly *_p = p_Copy(self._poly, _ring)
2382
2383	if fixed is not None:
2384	for m,v in fixed.iteritems():
2385	if PY_TYPE_CHECK(m,int) or PY_TYPE_CHECK(m,Integer):
2386	mi = m+1
2387	elif PY_TYPE_CHECK(m,MPolynomial_libsingular) and <MPolynomialRing_libsingular>m.parent() is parent:
2388	for i from 0 < i <= _ring.N:
2389	if p_GetExp((<MPolynomial_libsingular>m)._poly, i, _ring) != 0:
2390	mi = i
2391	break
2392	if i > _ring.N:
2393	raise TypeError, "key does not match"
2394	else:
2395	raise TypeError, "keys do not match self's parent"
2396	_p = pSubst(_p, mi, (<MPolynomial_libsingular>parent._coerce_c(v))._poly)
2397
2398	gd = parent.gens_dict()
2399	for m,v in kw.iteritems():
2400	m = gd[m]
2401	for i from 0 < i <= _ring.N:
2402	if p_GetExp((<MPolynomial_libsingular>m)._poly, i, _ring) != 0:
2403	mi = i
2404	break
2405	if i > _ring.N:
2406	raise TypeError, "key does not match"
2407	_p = pSubst(_p, mi, (<MPolynomial_libsingular>parent._coerce_c(v))._poly)
2408
2409	return new_MP(parent,_p)
2410
2411	def monomials(self):
2412	"""
2413	Return the list of monomials in self. The returned list is
2414	ordered by the term ordering of self.parent().
2415
2416	EXAMPLE:
2417	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
2418	sage: f = x + 3/2yz^2 + 2/3
2419	sage: f.monomials()
2420	[y*z^2, x, 1]
2421	sage: f = P(3/2)
2422	sage: f.monomials()
2423	[1]
2424
2425	TESTS:
2426	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
2427	sage: f = x
2428	sage: f.monomials()
2429	[x]
2430	sage: f = P(0)
2431	sage: f.monomials()
2432	[0]
2433
2434
2435	"""
2436	l = list()
2437	cdef MPolynomialRing_libsingular parent = <MPolynomialRing_libsingular>self._parent
2438	cdef ring *_ring = parent._ring
2439	cdef poly *p = p_Copy(self._poly, _ring)
2440	cdef poly *t
2441
2442	if p == NULL:
2443	return [parent._zero]
2444
2445	while p:
2446	t = pNext(p)
2447	p.next = NULL
2448	p_SetCoeff(p, n_Init(1,_ring), _ring)
2449	p_Setm(p, _ring)
2450	l.append( new_MP(parent,p) )
2451	p = t
2452
2453	return l
2454
2455	def constant_coefficient(self):
2456	"""
2457	Return the constant coefficient of this multivariate polynomial.
2458
2459	EXAMPLES:
2460	sage: P.<x, y> = MPolynomialRing(QQ,2)
2461	sage: f = 3x^2 - 2y + 7x^2y^2 + 5
2462	sage: f.constant_coefficient()
2463	5
2464	sage: f = 3*x^2
2465	sage: f.constant_coefficient()
2466	0
2467	"""
2468
2469	cdef poly *p = self._poly
2470	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2471	if p == NULL:
2472	return (<MPolynomialRing_libsingular>self._parent)._base._zero_element
2473
2474	while p.next:
2475	p = pNext(p)
2476
2477	if p_LmIsConstant(p, r):
2478	return co.si2sa( p_GetCoeff(p, r), r, (<MPolynomialRing_libsingular>self._parent)._base )
2479	else:
2480	return (<MPolynomialRing_libsingular>self._parent)._base._zero_element
2481
2482	def is_univariate(self):
2483	"""
2484	Return True if self is a univariate polynomial, that is if
2485	self contains only one variable.
2486
2487	EXAMPLE:
2488	sage: P.<x,y,z> = MPolynomialRing(GF(2),3)
2489	sage: f = x^2 + 1
2490	sage: f.is_univariate()
2491	True
2492	sage: f = y*x^2 + 1
2493	sage: f.is_univariate()
2494	False
2495	sage: f = P(0)
2496	sage: f.is_univariate()
2497	True
2498	"""
2499	return bool(len(self._variable_indices_(sort=False))<2)
2500
2501	def univariate_polynomial(self, R=None):
2502	"""
2503	Returns a univariate polynomial associated to this
2504	multivariate polynomial.
2505
2506	INPUT:
2507	R -- (default: None) PolynomialRing
2508
2509	If this polynomial is not in at most one variable, then a
2510	ValueError exception is raised. This is checked using the
2511	is_univariate() method. The new Polynomial is over the same
2512	base ring as the given MPolynomial and in the variable 'x' if
2513	no ring 'ring' is provided.
2514
2515	EXAMPLES:
2516	sage: R.<x, y> = MPolynomialRing(QQ,2)
2517	sage: f = 3x^2 - 2y + 7x^2y^2 + 5
2518	sage: f.univariate_polynomial()
2519	Traceback (most recent call last):
2520	...
2521	TypeError: polynomial must involve at most one variable
2522	sage: g = f.subs({x:10}); g
2523	700y^2 - 2y + 305
2524	sage: g.univariate_polynomial ()
2525	700x^2 - 2x + 305
2526	sage: g.univariate_polynomial(PolynomialRing(QQ,'z'))
2527	700z^2 - 2z + 305
2528	"""
2529	cdef poly *p = self._poly
2530	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2531	k = self.base_ring()
2532
2533	if not self.is_univariate():
2534	raise TypeError, "polynomial must involve at most one variable"
2535
2536	#construct ring if none
2537	if R == None:
2538	R = PolynomialRing(k,'x')
2539
2540	zero = k(0)
2541	coefficients = [zero] * (self.degree() + 1)
2542
2543	while p:
2544	coefficients[pTotaldegree(p, r)] = co.si2sa(p_GetCoeff(p, r), r, k)
2545	p = pNext(p)
2546
2547	return R(coefficients)
2548
2549
2550	def _variable_indices_(self, sort=True):
2551	"""
2552	Return the indices of all variables occuring in self.
2553	This index is the index as SAGE uses them (starting at zero), not
2554	as SINGULAR uses them (starting at one).
2555
2556	INPUT:
2557	sort -- specifies whether the indices shall be sorted
2558
2559	EXAMPLE:
2560	sage: P.<x,y,z> = MPolynomialRing(GF(2),3)
2561	sage: f = x*z^2 + z + 1
2562	sage: f._variable_indices_()
2563	[0, 2]
2564
2565	"""
2566	cdef poly *p
2567	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2568	cdef int i
2569	s = set()
2570	p = self._poly
2571	while p:
2572	for i from 1 <= i <= r.N:
2573	if p_GetExp(p,i,r):
2574	s.add(i-1)
2575	p = pNext(p)
2576	if sort:
2577	return sorted(s)
2578	else:
2579	return list(s)
2580
2581	def variables(self, sort=True):
2582	"""
2583	Return a list of all variables occuring in self.
2584
2585	INPUT:
2586	sort -- specifies whether the indices shall be sorted
2587
2588	EXAMPLE:
2589	sage: P.<x,y,z> = MPolynomialRing(GF(2),3)
2590	sage: f = x*z^2 + z + 1
2591	sage: f.variables()
2592	[z, x]
2593	sage: f.variables(sort=False)
2594	[x, z]
2595
2596	"""
2597	cdef poly p, v
2598	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2599	cdef int i
2600	l = list()
2601	si = set()
2602	p = self._poly
2603	while p:
2604	for i from 1 <= i <= r.N:
2605	if i not in si and p_GetExp(p,i,r):
2606	v = p_ISet(1,r)
2607	p_SetExp(v, i, 1, r)
2608	p_Setm(v, r)
2609	l.append(new_MP(self._parent, v))
2610	si.add(i)
2611	p = pNext(p)
2612	if sort:
2613	return sorted(l)
2614	else:
2615	return l
2616
2617	def variable(self, i=0):
2618	"""
2619	Return the i-th variable occuring in self. The index i is the
2620	index in self.variables().
2621
2622	EXAMPLE:
2623	sage: P.<x,y,z> = MPolynomialRing(GF(2),3)
2624	sage: f = x*z^2 + z + 1
2625	sage: f.variables()
2626	[z, x]
2627	sage: f.variable(1)
2628	x
2629	"""
2630	return self.variables()[i]
2631
2632	def nvariables(self):
2633	"""
2634	"""
2635	return self._variable_indices_(sort=False)
2636
2637	def is_constant(self):
2638	"""
2639	"""
2640	return bool(p_IsConstant(self._poly, (<MPolynomialRing_libsingular>self._parent)._ring))
2641
2642	cdef int is_constant_c(self):
2643	return p_IsConstant(self._poly, (<MPolynomialRing_libsingular>self._parent)._ring)
2644
2645	def __hash__(self):
2646	"""
2647	"""
2648	s = p_String(self._poly, (<MPolynomialRing_libsingular>self._parent)._ring, (<MPolynomialRing_libsingular>self._parent)._ring)
2649	return hash(s)
2650
2651	def lm(MPolynomial_libsingular self):
2652	"""
2653	Returns the lead monomial of self with respect to the term
2654	order of self.parent(). In SAGE a monomial is a product of
2655	variables in some power without a coefficient.
2656
2657	EXAMPLES:
2658	sage: R.<x,y,z>=PolynomialRing(GF(7),3,order='lex')
2659	sage: R.<x,y,z>=MPolynomialRing(GF(7),3,order='lex')
2660	sage: f = x^1y^2 + y^3z^4
2661	sage: f.lm()
2662	x*y^2
2663	sage: f = x^3y^2z^4 + x^3y^2z^1
2664	sage: f.lm()
2665	x^3y^2z^4
2666
2667	sage: R.<x,y,z>=MPolynomialRing(QQ,3,order='deglex')
2668	sage: f = x^1y^2z^3 + x^3y^2z^0
2669	sage: f.lm()
2670	xy^2z^3
2671	sage: f = x^1y^2z^4 + x^1y^1z^5
2672	sage: f.lm()
2673	xy^2z^4
2674
2675	sage: R.<x,y,z>=PolynomialRing(GF(127),3,order='degrevlex')
2676	sage: f = x^1y^5z^2 + x^4y^1z^3
2677	sage: f.lm()
2678	xy^5z^2
2679	sage: f = x^4y^7z^1 + x^4y^2z^3
2680	sage: f.lm()
2681	x^4y^7z
2682
2683	"""
2684	cdef poly *_p
2685	cdef ring *_ring
2686	_ring = (<MPolynomialRing_libsingular>self._parent)._ring
2687	_p = p_Head(self._poly, _ring)
2688	p_SetCoeff(_p, n_Init(1,_ring), _ring)
2689	p_Setm(_p,_ring)
2690	return new_MP((<MPolynomialRing_libsingular>self._parent), _p)
2691
2692
2693	def lc(MPolynomial_libsingular self):
2694	"""
2695	Leading coefficient of self. See self.lm() for details.
2696	"""
2697
2698	cdef poly *_p
2699	cdef ring *_ring
2700	cdef number *_n
2701	_ring = (<MPolynomialRing_libsingular>self._parent)._ring
2702
2703	if(_ring != currRing): rChangeCurrRing(_ring)
2704
2705	_p = p_Head(self._poly, _ring)
2706	_n = p_GetCoeff(_p, _ring)
2707
2708	return co.si2sa(_n, _ring, (<MPolynomialRing_libsingular>self._parent)._base)
2709
2710	def lt(MPolynomial_libsingular self):
2711	"""
2712	Leading term of self. In SAGE a term is a product of variables
2713	in some power AND a coefficient.
2714
2715	See self.lm() for details
2716	"""
2717	return new_MP((<MPolynomialRing_libsingular>self._parent),
2718	p_Head(self._poly,(<MPolynomialRing_libsingular>self._parent)._ring))
2719
2720	def is_zero(self):
2721	if self._poly is NULL:
2722	return True
2723	else:
2724	return False
2725
2726	def __nonzero__(self):
2727	if self._poly:
2728	return True
2729	else:
2730	return False
2731
2732	def __floordiv__(self, right):
2733	"""
2734	Perform division with remainder and return the quotient.
2735
2736	INPUT:
2737	right -- something coercable to an MPolynomial_libsingular in self.parent()
2738
2739	EXAMPLES:
2740	sage: R.<x,y,z> = PolynomialRing(GF(32003),3)
2741	sage: R.<x,y,z> = MPolynomialRing(GF(32003),3)
2742	sage: f = y*x^2 + x + 1
2743	sage: f//x
2744	x*y + 1
2745	sage: f//y
2746	x^2
2747	"""
2748	cdef MPolynomialRing_libsingular parent = <MPolynomialRing_libsingular>(<MPolynomial_libsingular>self)._parent
2749	cdef ring *r = parent._ring
2750	if(r != currRing): rChangeCurrRing(r)
2751	cdef MPolynomial_libsingular _self, _right
2752	cdef poly *quo
2753
2754	_self = self
2755
2756	if not PY_TYPE_CHECK(right, MPolynomial_libsingular) or (<ParentWithBase>parent is not (<MPolynomial_libsingular>right)._parent):
2757	_right = parent._coerce_c(right)
2758	else:
2759	_right = right
2760
2761	if right.is_zero():
2762	raise ZeroDivisionError
2763
2764	quo = singclap_pdivide( _self._poly, _right._poly )
2765	return new_MP(parent, quo)
2766
2767	def factor(self, param=0):
2768	"""
2769	Return the factorization of self.
2770
2771	INPUT:
2772	param -- 0: returns factors and multiplicities, first factor is a constant.
2773	1: returns non-constant factors (no multiplicities).
2774	2: returns non-constant factors and multiplicities.
2775	EXAMPLE:
2776	sage: R.<x,y,z> = PolynomialRing(GF(32003),3)
2777	sage: R.<x,y,z> = MPolynomialRing(GF(32003),3)
2778	sage: f = 9(x-1)^2(y+z)
2779	sage: f.factor(0)
2780	9 * (y + z) * (x - 1)^2
2781	sage: f.factor(1)
2782	(y + z) * (x - 1)
2783	sage: f.factor(2)
2784	(y + z) * (x - 1)^2
2785
2786	"""
2787	cdef ring *_ring
2788	cdef intvec *iv
2789	cdef int *ivv
2790	cdef ideal *I
2791	cdef MPolynomialRing_libsingular parent
2792	cdef int i
2793
2794	parent = self._parent
2795	_ring = parent._ring
2796
2797	if(_ring != currRing): rChangeCurrRing(_ring)
2798
2799	iv = NULL
2800	I = singclap_factorize ( self._poly, &iv , int(param)) #delete iv at some point
2801
2802	if param==1:
2803	v = [(new_MP(parent, p_Copy(I.m[i],_ring)) , 1) for i in range(I.ncols)]
2804	else:
2805	ivv = iv.ivGetVec()
2806	v = [(new_MP(parent, p_Copy(I.m[i],_ring)) , ivv[i]) for i in range(I.ncols)]
2807	oo = (new_MP(parent, p_ISet(1,_ring)),1)
2808	if oo in v:
2809	v.remove(oo)
2810
2811	F = Factorization(v)
2812	F.sort()
2813
2814	delete(iv)
2815	id_Delete(&I,_ring)
2816
2817	return F
2818
2819	def lift(self, I):
2820	"""
2821	given an ideal I = (f_1,...,f_r) and some g (== self) in I,
2822	find s_1,...,s_r such that g = s_1 f_1 + ... + s_r f_r
2823
2824	EXAMPLE:
2825	sage: A.<x,y> = PolynomialRing(QQ,2,order='degrevlex')
2826	sage: A.<x,y> = MPolynomialRing(QQ,2,order='degrevlex')
2827	sage: I = A.ideal([x^10 + x^9y^2, y^8 - x^2y^7 ])
2828	sage: f = x*y^13 + y^12
2829	sage: M = f.lift(I)
2830	sage: M
2831	[y^7, x^7y^2 + x^8 + x^5y^3 + x^6y + x^3y^4 + x^4y^2 + xy^5 + x^2*y^3 + y^4]
2832	sage: sum( map( mul , zip( M, I.gens() ) ) ) == f
2833	True
2834	"""
2835
2836	cdef ideal *fI = idInit(1,1)
2837	cdef ideal *_I
2838	cdef MPolynomialRing_libsingular parent = <MPolynomialRing_libsingular>self._parent
2839	cdef int i = 0
2840	cdef int j
2841	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2842	cdef ideal *res
2843
2844	if PY_TYPE_CHECK(I, MPolynomialIdeal):
2845	I = I.gens()
2846
2847	_I = idInit(len(I),1)
2848
2849	for f in I:
2850	if not (PY_TYPE_CHECK(f,MPolynomial_libsingular) \
2851	and <MPolynomialRing_libsingular>(<MPolynomial_libsingular>f)._parent is parent):
2852	try:
2853	f = parent._coerce_c(f)
2854	except TypeError, msg:
2855	id_Delete(&fI,r)
2856	id_Delete(&_I,r)
2857	raise TypeError, msg
2858
2859	_I.m[i] = p_Copy((<MPolynomial_libsingular>f)._poly, r)
2860	i+=1
2861
2862	fI.m[0]= p_Copy(self._poly, r)
2863
2864	res = idLift(_I, fI, NULL, 0, 0, 0)
2865	l = []
2866	for i from 0 <= i < IDELEMS(res):
2867	for j from 1 <= j <= IDELEMS(_I):
2868	l.append( new_MP(parent, pTakeOutComp1(&res.m[i], j)) )
2869
2870	id_Delete(&fI, r)
2871	id_Delete(&_I, r)
2872	id_Delete(&res, r)
2873	return Sequence(l, check=False, immutable=True)
2874
2875	def reduce(self,I):
2876	"""
2877	Return the normal form of self w.r.t. I, i.e. return the
2878	remainder of self with respect to the polynomials in I. If the
2879	polynomial set/list I is not a Groebner basis the result is
2880	not canonical.
2881
2882	INPUT:
2883	I -- a list/set of polynomials in self.parent(). If I is an ideal,
2884	the generators are used.
2885
2886	EXAMPLE:
2887
2888	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
2889	sage: f1 = -2 * x^2 + x^3
2890	sage: f2 = -2 * y + x* y
2891	sage: f3 = -x^2 + y^2
2892	sage: F = Ideal([f1,f2,f3])
2893	sage: g = xy - 3x*y^2
2894	sage: g.reduce(F)
2895	-6y^2 + 2y
2896	sage: g.reduce(F.gens())
2897	-6y^2 + 2y
2898
2899	"""
2900	cdef ideal *_I
2901	cdef MPolynomialRing_libsingular parent = <MPolynomialRing_libsingular>self._parent
2902	cdef int i = 0
2903	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
2904	cdef poly *res
2905
2906	if(r != currRing): rChangeCurrRing(r)
2907
2908	if PY_TYPE_CHECK(I, MPolynomialIdeal):
2909	I = I.gens()
2910
2911	_I = idInit(len(I),1)
2912	for f in I:
2913	if not (PY_TYPE_CHECK(f,MPolynomial_libsingular) \
2914	and <MPolynomialRing_libsingular>(<MPolynomial_libsingular>f)._parent is parent):
2915	try:
2916	f = parent._coerce_c(f)
2917	except TypeError, msg:
2918	id_Delete(&_I,r)
2919	raise TypeError, msg
2920
2921	_I.m[i] = p_Copy((<MPolynomial_libsingular>f)._poly, r)
2922	i+=1
2923
2924	#the second parameter would be qring!
2925	res = kNF(_I, NULL, p_Copy(self._poly, r))
2926	return new_MP(parent,res)
2927
2928	def gcd(self, right):
2929	"""
2930	Return the greates common divisor of self and right.
2931
2932	INPUT:
2933	right -- polynomial
2934
2935	EXAMPLES:
2936	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
2937	sage: f = (xyz)^6 - 1
2938	sage: g = (xyz)^4 - 1
2939	sage: f.gcd(g)
2940	x^2y^2z^2 - 1
2941	sage: GCD([x^3 - 3*x + 2, x^4 - 1, x^6 -1])
2942	x - 1
2943
2944	TESTS:
2945	sage: Q.<x,y,z> = MPolynomialRing(QQ,3)
2946	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
2947	sage: P(0).gcd(Q(0))
2948	0
2949	sage: x.gcd(1)
2950	1
2951
2952	"""
2953	cdef MPolynomial_libsingular _right
2954	cdef poly *_res
2955	cdef ring *_ring
2956
2957	_ring = (<MPolynomialRing_libsingular>self._parent)._ring
2958
2959	if(_ring != currRing): rChangeCurrRing(_ring)
2960
2961	if not PY_TYPE_CHECK(right, MPolynomial_libsingular):
2962	_right = (<MPolynomialRing_libsingular>self._parent)._coerce_c(right)
2963	else:
2964	_right = (<MPolynomial_libsingular>right)
2965
2966	_res = singclap_gcd(p_Copy(self._poly, _ring), p_Copy(_right._poly, _ring))
2967
2968	return new_MP((<MPolynomialRing_libsingular>self._parent), _res)
2969
2970	def lcm(self, MPolynomial_libsingular g):
2971	"""
2972	Return the least common multiple of self and g.
2973
2974	INPUT:
2975	g -- polynomial
2976
2977	OUTPUT:
2978	polynomial
2979
2980	EXAMPLE:
2981	sage: P.<x,y,z> = MPolynomialRing(QQ,3)
2982	sage: p = (x+y)*(y+z)
2983	sage: q = (z^4+2)*(y+z)
2984	sage: lcm(p,q)
2985	xyz^4 + y^2z^4 + xz^5 + yz^5 + 2xy + 2y^2 + 2xz + 2yz
2986
2987	NOTE: This only works for GF(p) and QQ as base rings
2988	"""
2989	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
2990	cdef poly ret, prod, *gcd
2991	if(_ring != currRing): rChangeCurrRing(_ring)
2992
2993	if self._parent is not g._parent:
2994	g = (<MPolynomialRing_libsingular>self._parent)._coerce_c(g)
2995
2996	gcd = singclap_gcd(p_Copy(self._poly, _ring), p_Copy((<MPolynomial_libsingular>g)._poly, _ring))
2997	prod = pp_Mult_qq(self._poly, (<MPolynomial_libsingular>g)._poly, _ring)
2998	ret = singclap_pdivide(prod , gcd )
2999	p_Delete(&prod, _ring)
3000	p_Delete(&gcd, _ring)
3001	return new_MP(self._parent, ret)
3002
3003	def is_square_free(self):
3004	"""
3005	"""
3006	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
3007	if(_ring != currRing): rChangeCurrRing(_ring)
3008	return bool(singclap_isSqrFree(self._poly))
3009
3010	def quo_rem(self, MPolynomial_libsingular right):
3011	"""
3012	Returns quotient and remainder of self and right.
3013
3014	EXAMPLES:
3015	sage: R.<x,y> = MPolynomialRing(QQ,2)
3016	sage: f = y*x^2 + x + 1
3017	sage: f.quo_rem(x)
3018	(x*y + 1, 1)
3019	sage: f.quo_rem(y)
3020	(x^2, x + 1)
3021
3022	TESTS:
3023	sage: R.<x,y> = MPolynomialRing(QQ,2)
3024	sage: R(0).quo_rem(R(1))
3025	(0, 0)
3026	sage: R(1).quo_rem(R(0))
3027	Traceback (most recent call last):
3028	...
3029	ZeroDivisionError
3030
3031	"""
3032	#cdef ideal selfI, rightI, R, res
3033	cdef poly quo, rem
3034	cdef MPolynomialRing_libsingular parent = <MPolynomialRing_libsingular>self._parent
3035	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
3036	if(r != currRing): rChangeCurrRing(r)
3037
3038	if self._parent is not right._parent:
3039	right = self._parent._coerce_c(right)
3040
3041	if right.is_zero():
3042	raise ZeroDivisionError
3043
3044	quo = singclap_pdivide( self._poly, right._poly )
3045	rem = p_Add_q(p_Copy(self._poly, r), p_Neg(pp_Mult_qq(right._poly, quo, r), r), r)
3046	return new_MP(parent, quo), new_MP(parent, rem)
3047
3048	def _magma_(self, magma=None):
3049	"""
3050	Returns the MAGMA representation of self.
3051
3052	EXAMPLES:
3053	sage: R.<x,y> = MPolynomialRing(GF(2),2)
3054	sage: f = y*x^2 + x +1
3055	sage: f._magma_() #optional
3056	x^2*y + x + 1
3057	"""
3058	if magma is None:
3059	# TODO: import this globally
3060	import sage.interfaces.magma
3061	magma = sage.interfaces.magma.magma
3062
3063	magma_gens = [e.name() for e in self.parent()._magma_().gens()]
3064	f = self._repr_with_changed_varnames(magma_gens)
3065	return magma(f)
3066
3067	def _singular_(self, singular=singular_default, have_ring=False):
3068	"""
3069	Return a SINGULAR (as in the CAS) element for this
3070	element. The result is cached.
3071
3072	INPUT:
3073	singular -- interpreter (default: singular_default)
3074	have_ring -- should the correct ring not be set in SINGULAR first (default:False)
3075
3076	EXAMPLES:
3077	sage: P.<x,y,z> = PolynomialRing(GF(127),3)
3078	sage: x._singular_()
3079	x
3080	sage: f =(x^2 + 35*y + 128); f
3081	x^2 + 35*y + 1
3082	sage: x._singular_().name() == x._singular_().name()
3083	True
3084
3085
3086	TESTS:
3087	sage: P.<x,y,z> = MPolynomialRing(GF(127),3)
3088	sage: P.<x,y,z> = PolynomialRing(GF(127),3)
3089	sage: P(0)._singular_()
3090	0
3091
3092	"""
3093	if not have_ring:
3094	self.parent()._singular_(singular).set_ring() #this is expensive
3095
3096	try:
3097	if self.__singular is None:
3098	return self._singular_init_c(singular, True)
3099
3100	self.__singular._check_valid()
3101
3102	if self.__singular.parent() is singular:
3103	return self.__singular
3104
3105	except (AttributeError, ValueError):
3106	pass
3107
3108	return self._singular_init_c(singular, True)
3109
3110	def _singular_init_(self,singular=singular_default, have_ring=False):
3111	"""
3112	Return a new SINGULAR (as in the CAS) element for this element.
3113
3114	INPUT:
3115	singular -- interpreter (default: singular_default)
3116	have_ring -- should the correct ring not be set in SINGULAR first (default:False)
3117
3118	EXAMPLES:
3119	sage: P.<x,y,z> = PolynomialRing(GF(127),3)
3120	sage: x._singular_init_()
3121	x
3122	sage: (x^2+37*y+128)._singular_init_()
3123	x^2+37*y+1
3124	sage: x._singular_init_().name() == x._singular_init_().name()
3125	False
3126
3127	TESTS:
3128	sage: P(0)._singular_init_()
3129	0
3130	"""
3131	return self._singular_init_c(singular, have_ring)
3132
3133	cdef _singular_init_c(self,singular, have_ring):
3134	"""
3135	See MPolynomial_libsingular._singular_init_
3136
3137	"""
3138	if not have_ring:
3139	self.parent()._singular_(singular).set_ring() #this is expensive
3140
3141	self.__singular = singular(str(self))
3142	return self.__singular
3143
3144	def sub_m_mul_q(self, MPolynomial_libsingular m, MPolynomial_libsingular q):
3145	"""
3146	Return self - m*q, where m must be a monomial and q a
3147	polynomial.
3148
3149	INPUT:
3150	m -- a monomial
3151	q -- a polynomial
3152
3153	EXAMPLE:
3154	sage: P.<x,y,z>=PolynomialRing(QQ,3)
3155	sage: x.sub_m_mul_q(y,z)
3156	-y*z + x
3157
3158	TESTS:
3159	sage: from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
3160	sage: Q.<x,y,z>=MPolynomialRing(QQ,3)
3161	sage: P.<x,y,z>=MPolynomialRing(QQ,3)
3162	sage: P(0).sub_m_mul_q(P(0),P(1))
3163	0
3164	sage: x.sub_m_mul_q(Q.gen(1),Q.gen(2))
3165	-y*z + x
3166
3167	"""
3168	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
3169
3170	if not self._parent is m._parent:
3171	m = self._parent._coerce_c(m)
3172	if not self._parent is q._parent:
3173	q = self._parent._coerce_c(q)
3174
3175	if m._poly and m._poly.next:
3176	raise ArithmeticError, "m must be a monomial"
3177	elif not m._poly:
3178	return self
3179
3180	return new_MP(self._parent, p_Minus_mm_Mult_qq(p_Copy(self._poly, r), m._poly, q._poly, r))
3181
3182	def add_m_mul_q(self, MPolynomial_libsingular m, MPolynomial_libsingular q):
3183	"""
3184	Return self + m*q, where m must be a monomial and q a
3185	polynomial.
3186
3187	INPUT:
3188	m -- a monomial
3189	q -- a polynomial
3190
3191	EXAMPLE:
3192	sage: P.<x,y,z>=PolynomialRing(QQ,3)
3193	sage: x.add_m_mul_q(y,z)
3194	y*z + x
3195
3196	TESTS:
3197	sage: from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
3198	sage: R.<x,y,z>=MPolynomialRing(QQ,3)
3199	sage: P.<x,y,z>=MPolynomialRing(QQ,3)
3200	sage: P(0).add_m_mul_q(P(0),P(1))
3201	0
3202	sage: x.add_m_mul_q(R.gen(),R.gen(1))
3203	x*y + x
3204	"""
3205
3206	cdef ring *r = (<MPolynomialRing_libsingular>self._parent)._ring
3207
3208	if not self._parent is m._parent:
3209	m = self._parent._coerce_c(m)
3210	if not self._parent is q._parent:
3211	q = self._parent._coerce_c(q)
3212
3213	if m._poly and m._poly.next:
3214	raise ArithmeticError, "m must be a monomial"
3215	elif not m._poly:
3216	return self
3217
3218	return new_MP(self._parent, p_Plus_mm_Mult_qq(p_Copy(self._poly, r), m._poly, q._poly, r))
3219
3220
3221	def __reduce__(self):
3222	"""
3223
3224	Serialize self.
3225
3226	EXAMPLES:
3227	sage: P.<x,y,z> = PolynomialRing(QQ,3, order='degrevlex')
3228	sage: f = 27/113 * x^2 + y*z + 1/2
3229	sage: f == loads(dumps(f))
3230	True
3231
3232	sage: P = PolynomialRing(GF(127),3,names='abc')
3233	sage: a,b,c = P.gens()
3234	sage: f = 57 * a^2b + 43 c + 1
3235	sage: f == loads(dumps(f))
3236	True
3237
3238	"""
3239	return sage.rings.polynomial.multi_polynomial_libsingular.unpickle_MPolynomial_libsingular, ( self._parent, self.dict() )
3240
3241	def _im_gens_(self, codomain, im_gens):
3242	"""
3243
3244	INPUT:
3245	codomain
3246	im_gens
3247
3248	EXAMPLES:
3249	sage: R.<x,y> = PolynomialRing(QQ, 2)
3250	sage: f = R.hom([y,x], R)
3251	sage: f(x^2 + 3*y^5)
3252	3*x^5 + y^2
3253	"""
3254	#TODO: very slow
3255	n = self.parent().ngens()
3256	if n == 0:
3257	return codomain._coerce_(self)
3258	y = codomain(0)
3259	for (m,c) in self.dict().iteritems():
3260	y += codomain(c)mul([ im_gens[i]*m[i] for i in range(n) ])
3261	return y
3262
3263	def diff(self, MPolynomial_libsingular variable, have_ring=True):
3264	"""
3265	Differentiates self with respect to the provided variable. This
3266	is completely symbolic so it is also defined over e.g. finite
3267	fields.
3268
3269	INPUT:
3270	variable -- the derivative is taken with respect to variable
3271	have_ring -- ignored, accepted for compatibility reasons
3272
3273	EXAMPLES:
3274	sage: R.<x,y> = PolynomialRing(QQ,2)
3275	sage: f = 3x^3y^2 + 5y^2 + 3x + 2
3276	sage: f.diff(x)
3277	9x^2y^2 + 3
3278	sage: f.diff(y)
3279	6x^3y + 10*y
3280
3281	The derivate is also defined over finite fields:
3282
3283	sage: R.<x,y> = PolynomialRing(GF(2**8, 'a'),2)
3284	sage: f = x^3*y^2 + y^2 + x + 2
3285	sage: f.diff(x)
3286	x^2*y^2 + 1
3287
3288	"""
3289	cdef int i, var_i
3290
3291	cdef poly *p
3292	if variable._parent is not self._parent:
3293	raise TypeError, "provided variable is not in same ring as self"
3294	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
3295	if _ring != currRing: rChangeCurrRing(_ring)
3296
3297	var_i = -1
3298	for i from 0 <= i <= _ring.N:
3299	if p_GetExp(variable._poly, i, _ring):
3300	if var_i == -1:
3301	var_i = i
3302	else:
3303	raise TypeError, "provided variable is not univariate"
3304
3305	if var_i == -1:
3306	raise TypeError, "provided variable is constant"
3307
3308
3309	p = pDiff(self._poly, var_i)
3310	return new_MP(self._parent,p)
3311
3312	def resultant(self, MPolynomial_libsingular other, variable=None):
3313	"""
3314	computes the resultant of self and the first argument with
3315	respect to the variable given as the second argument.
3316
3317	If a second argument is not provide the first variable of
3318	self.parent() is chosen.
3319
3320	INPUT:
3321	other -- polynomial in self.parent()
3322	variable -- optional variable (of type polynomial) in self.parent() (default: None)
3323
3324	EXAMPLE:
3325	sage: P.<x,y> = PolynomialRing(QQ,2)
3326	sage: a = x+y
3327	sage: b = x^3-y^3
3328	sage: c = a.resultant(b); c
3329	-2*y^3
3330	sage: d = a.resultant(b,y); d
3331	2*x^3
3332
3333	The SINGULAR example:
3334	sage: R.<x,y,z> = PolynomialRing(GF(32003),3)
3335	sage: f = 3 * (x+2)^3 + y
3336	sage: g = x+y+z
3337	sage: f.resultant(g,x)
3338	3y^3 + 9y^2z + 9yz^2 + 3z^3 - 18y^2 - 36yz - 18z^2 + 35y + 36z - 24
3339
3340	TESTS:
3341	sage: from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
3342	sage: P.<x,y> = MPolynomialRing_libsingular(QQ,2,order='degrevlex')
3343	sage: a = x+y
3344	sage: b = x^3-y^3
3345	sage: c = a.resultant(b); c
3346	-2*y^3
3347	sage: d = a.resultant(b,y); d
3348	2*x^3
3349
3350	"""
3351	cdef ring *_ring = (<MPolynomialRing_libsingular>self._parent)._ring
3352	cdef poly *rt
3353
3354	if variable is None:
3355	variable = self.parent().gen(0)
3356
3357	if not self._parent is other._parent:
3358	raise TypeError, "first parameter needs to be an element of self.parent()"
3359
3360	if not variable.parent() is self.parent():
3361	raise TypeError, "second parameter needs to be an element of self.parent() or None"
3362
3363	rt = singclap_resultant(self._poly, other._poly, (<MPolynomial_libsingular>variable)._poly )
3364	return new_MP(self._parent, rt)
3365
3366	def unpickle_MPolynomial_libsingular(MPolynomialRing_libsingular R, d):
3367	"""
3368	Deserialize a MPolynomial_libsingular object
3369
3370	INPUT:
3371	R -- the base ring
3372	d -- a Python dictionary as returned by MPolynomial_libsingular.dict
3373
3374	"""
3375	cdef ring *r = R._ring
3376	cdef poly m, p
3377	cdef int _i, _e
3378	p = p_ISet(0,r)
3379	rChangeCurrRing(r)
3380	for mon,c in d.iteritems():
3381	m = p_Init(r)
3382	for i,e in mon.sparse_iter():
3383	_i = i
3384	if _i >= r.N:
3385	p_Delete(&p,r)
3386	p_Delete(&m,r)
3387	raise TypeError, "variable index too big"
3388	_e = e
3389	if _e <= 0:
3390	p_Delete(&p,r)
3391	p_Delete(&m,r)
3392	raise TypeError, "exponent too small"
3393	p_SetExp(m, _i+1,_e, r)
3394	p_SetCoeff(m, co.sa2si(c, r), r)
3395	p_Setm(m,r)
3396	p = p_Add_q(p,m,r)
3397	return new_MP(R,p)
3398
3399
3400	cdef poly addwithcarry(poly tempvector, poly maxvector, int pos, ring _ring):
3401	if p_GetExp(tempvector, pos, _ring) < p_GetExp(maxvector, pos, _ring):
3402	p_SetExp(tempvector, pos, p_GetExp(tempvector, pos, _ring)+1, _ring)
3403	else:
3404	p_SetExp(tempvector, pos, 0, _ring)
3405	tempvector = addwithcarry(tempvector, maxvector, pos + 1, _ring)
3406	p_Setm(tempvector, _ring)
3407	return tempvector
3408
3409

Note: See TracBrowser for help on using the repository browser.

Download in other formats: