ASKSAGE: Sage Q&A Forum - Latest question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 06 Feb 2019 00:49:28 -0600Groebner basis computationhttp://ask.sagemath.org/question/45340/groebner-basis-computation/ How we can use SageMath for computing Groebner basis? e.g. we have 35 multivariate polynomial equations, is it possible to compute groebner basis at a time?MMumtazWed, 06 Feb 2019 00:49:28 -0600http://ask.sagemath.org/question/45340/Multiplying Roots of a Polynomialhttp://ask.sagemath.org/question/45077/multiplying-roots-of-a-polynomial/ Hello - I'm new to SAGE, so trying to get to grips with the basics!
I have a polynomial $f = x^{2}((x+\frac{1}{x})^{2} - 1) $. I would like to be able to do to things:
1) See $f$ in the form $\beta(x - \alpha_1)(x-\alpha_2)(x-\alpha_3)(x-\alpha_4)$
2) Know the value of $\beta (\prod_i \alpha_i)$
Using
>f.factor
>f.roots
Does not give me the linear factorisation, nor can I see a simple way of getting the product of the roots.
In the long term I'm hoping to do this for polynomials of larger degree, so any advice would be great!LukeC93Thu, 17 Jan 2019 05:29:29 -0600http://ask.sagemath.org/question/45077/resolution stuck system polynomials order twohttp://ask.sagemath.org/question/44663/resolution-stuck-system-polynomials-order-two/I have three variables (ISP, ch, $n_C$). I have two equations. My aim is to solve a system of two equations to compute the values of ISP and ch according to $n_C$.
My two equations are polynomials of order two. I put the equations in a framabin :
https://framabin.org/p/?2ea1412d13d7f562#vh/+dEhivcOqEQAl3RerwPLb4XPAYRB4ehQ+eAL65MM=
I tried the solvers of sage and simpy but it does not work. I tried to use Groebner basis but I failed.borostackThu, 13 Dec 2018 14:47:04 -0600http://ask.sagemath.org/question/44663/Conversion from symbolic expression to polynomial stuckhttp://ask.sagemath.org/question/44302/conversion-from-symbolic-expression-to-polynomial-stuck/I have a very long symbolic expression with six variables:
> (E_mu, E_xi3, ISP, T, V_mu, V_xi3, Z,
> m, mu, n_I, n_P, xi_1, xi_3)
I call the expression AN. I want to convert it in a polynomial of two variables (mu and xi_3). I tried the following command:
AP = AN.polynomial(None,ring=SR['mu,xi_3'])
This works for simple expressions but the command remains stuck for the given expression. My aim is to get the monomials of the polynomial. mu and xi_3 are random variables and I want to compute the expectation of AN (E[AN]). Thus, I could substitute the expectations in the expression.
AN is a numerator, so it is not a fraction. Moreover, I also tried simply_rational but it did not help to convert.
Is there anyway to convert easily ?
EDIT : [I put the symbolic expression in a framabin.](https://framabin.org/p/?2417942055fc9fc9#Oenpb+BI567QXYcVCzjO4zkLRC0byQaNdhi2XW/CUkY=)borostackFri, 16 Nov 2018 05:53:39 -0600http://ask.sagemath.org/question/44302/Implementing the AKS primality testhttp://ask.sagemath.org/question/44275/implementing-the-aks-primality-test/I am implementing the infamous AKS deterministic primality test using SageMath.
My code is as follows:
# AKS Primality Test
def fast_exponentation(a, n):
"""
Computes a^n fast
"""
ans = 1
while n:
if n & 1 : ans = ans * a
a = a * a
n >>= 1
return ans
# Determines whether n is a power a ^ b, b > 1
def is_perfect_power(n):
lgn = 1 + ( len( bin ( abs ( n ) ) ) - 2)
for b in range(2,lgn):
# we use binary search to check if n root b is an integer
lowa = 0
higha = n
while lowa < higha - 1:
mida = (lowa + higha) // 2
ab = fast_exponentation(mida,b)
if ab > n: higha = mida
elif ab < n: lowa = mida
else: return True # mida ^ b
return False
def AKS_primality_test(n):
if n == 1: return None
if n == 2: return True
# if n is of the form a^b for a > 1 and b > 1 return False
if is_perfect_power(n): return False
print("Passed perfect power test")
# find the smallest integer r such that
# a. gcd(n,r) > 1 or
# b. [n]_r has multiplicative order > 4 log2(n)^2
for r in xrange(2,n):
# if r shares a non trivial factor with n, n is composite
if gcd(r,n) > 1: return False
# compute the multiplicative order of n in Z mod r
ord_r = mod(n,r).multiplicative_order()
if ord_r > 4*log(n,2)^2: break
print("min r found! r={}".format(r))
if r == n: return True
# We get the ring where we will check for primality
ZnX = IntegerModRing(n)[x]
ZnX_div_f = ZnX.quo(x^r-1)
for j in range(1, 2*log(n,2)*sqrt(r)+2):
if ZnX_div_f((x+j)^n) != ZnX_div_f(x^n+j):
print("The freshman dream test failed for j={}".format(j))
return False
print("Passed final test")
return True
print("271 is prime = " + str(AKS_primality_test(271)))
For the implementation I am following the scheme in page 550 of *A Computational Introduction to Number Theory
and Algebra*, by Victor Shoup ([available online in his webpage](https://shoup.net/ntb/ntb-v2.pdf)).
![image description](/upfiles/1542188648972937.png)
However, the code above, which should test whether 271 is prime, gives us the following output:
Passed perfect power test
min r found! r=269
The freshman dream test failed for j=1
271 is prime = False
This is wrong - 271 is in fact prime. So I assume I am misunderstanding how to use quotient rings in SageMath.
Why does the equality test fail?
EDIT: as pointed out by @rburing, the line `if ZnX_div_f((x+j)^n) != ZnX_div_f(x+j^n): ` should read `if ZnX_div_f((x+j)^n) != ZnX_div_f(x^n+j): ` I have edited the code above to reflect this, but the behavior of the program is the same nonetheless.JsevillamolWed, 14 Nov 2018 03:50:59 -0600http://ask.sagemath.org/question/44275/How to define this polynomial?http://ask.sagemath.org/question/44061/how-to-define-this-polynomial/$ h_i(x_1,\ldots,x_n) = \sum_{j_1+\cdots+j_n=i} \;\prod_{k=1}^n x_k^{j_k}$
I think I need a Partitions(i) to do it but I need some $j_k$ to be zero, then I tried PartitionTuples(level=n,size=i), and it still doesn't work.GoGoKoWed, 24 Oct 2018 15:48:34 -0500http://ask.sagemath.org/question/44061/polynomial list, arrayhttp://ask.sagemath.org/question/11065/polynomial-list-array/I would like to store polynomials in an array/list in Sage. How do I do this?
-------
I am editing my question due to the confusion expressed below.
I would like to generate a list of polynomials from a loop. HansThu, 19 Jun 2014 15:30:28 -0500http://ask.sagemath.org/question/11065/Show correct output of polynomialhttp://ask.sagemath.org/question/43792/show-correct-output-of-polynomial/ Hi, I am new to sage and I am trying to construct a polynomial to try to transform it. However, when writing it I am getting an incorrect output.
x=var("x")
s = (x^2+2*x+1) + 1/(x^2+2*x+1)
s.show()
And this is the output I am getting -
x^2+2x + 1/x^2+2x+1 +1 (I am unsure why the 1 is carrying over all the way to right side)
I am looking to obtain the following output to begin transforming it -
x^2+2x+1 + 1/x^2+2x+1
Appreciate any help!nevar123Sat, 29 Sep 2018 14:15:10 -0500http://ask.sagemath.org/question/43792/defining boolean variables in sagehttp://ask.sagemath.org/question/26662/defining-boolean-variables-in-sage/ Hi Guys,
By writing this:
B.<a,b> = BooleanPolynomialRing()
Not only a Boolean Polynomial Ring in 'a' and 'b' is defined, but 'a' and 'b' are also treated as boolean variables.
However, if we write in this manner:
B = BooleanPolynomialRing(names = ['a','b'])
We'll obtain a Boolean Polynomial Ring in 'a' and 'b', but we don't even get 'a' and 'b' as variables.
Is there any way to resolve the issue in the second method, especially, if we have a boolean polynomial ring of >1000 variables? Thanks in advance!freako89Tue, 28 Apr 2015 01:56:42 -0500http://ask.sagemath.org/question/26662/Imaginary part of a ratio of polynomialshttp://ask.sagemath.org/question/42609/imaginary-part-of-a-ratio-of-polynomials/When working with polynomial rings with complex coefficients, how to pass the assumption that the unknown variable belongs to the Reals set?
I am trying to recover and plot the imaginary part of the the ratio of polynomials *arg2*, but as sage assumes the independent variable is complex, it takes a long time and returns a huge answer. But I think the answer could be much simpler if it assumes that the independent variable belongs to the Reals set. The *assume()* function doesn't seems to work in this case.
*arg* has complex terms in the numerator and denominator. *arg2* has complex terms only on the numerator. So, it should be easy to obtain the imaginary part of *arg2*.
s = polygen(CC, "s")
omega = polygen(CC, "omega")
I = CC(I)
# assume(omega, "real") # this does not work
s21 = 0.894292331092066/(s^4 + 2.14081977623463*s^3 + 3.15237117897600*s^2 + 2.31898630138664*s + 0.902488008823108)
s21_omega = s21.subs(s=I*omega)
s21_omega_diff = s21_omega.derivative()
arg = 1/s21_omega*s21_omega_diff
num = arg.numerator()
den = arg.denominator()
den_conj = den.map_coefficients(lambda z: z.conjugate())
num2 = (num*den_conj).map_coefficients(lambda z: 0 if abs(z.real())<1e-10 and abs(z.imag())<1e-10 else z.imag()*I if abs(z.real())<1e-10 and abs(z.imag())>=1e-10 else z.real() if abs(z.real())>=1e-10 and abs(z.imag())<1e-10 else z)
den2 = (den*den_conj).map_coefficients(lambda z: 0 if abs(z.real())<1e-10 and abs(z.imag())<1e-10 else z.imag()*I if abs(z.real())<1e-10 and abs(z.imag())>=1e-10 else z.real() if abs(z.real())>=1e-10 and abs(z.imag())<1e-10 else z)
arg2 = num2/den2
print(arg2)
print(imag(arg2))
# plot(-imag(arg2), (-1.5, 1.5)).show()
Output:
(-3.19903509380033*omega^15 - 1.71213939841908*I*omega^14 + 9.63823791915898*omega^13 + 3.11426578979509*I*omega^12 - 15.8129908994689*omega^11 - 4.60245136409720*I*omega^10 + 13.7326241083798*omega^9 + 1.24770230131999*I*omega^8 - 9.58497293293397*omega^7 - 0.760732566992820*I*omega^6 + 4.72291963120453*omega^5 - 2.52135706735033*I*omega^4 - 2.44038907746050*omega^3 - 0.463630243025939*I*omega^2 + 0.203401406783201*omega - 1.36326886734867*I)/(0.799758773450081*omega^16 - 2.75378226261685*omega^14 + 5.27099696648963*omega^12 - 5.49304964335183*omega^10 + 4.79248646646703*omega^8 - 3.14861308746971*omega^6 + 2.44038907746050*omega^4 - 0.406802813566401*omega^2 + 0.530548112702672)
2.55845638284151*imag_part(omega)^31/(0.639614095710378*imag_part(omega)^32 + 10.2338255313661*imag_part(omega)^30*real_part(omega)^2 + 76.7536914852462*imag_part(omega)^28*real_part(omega)^4 + 358.183893597743*imag_part(omega)^26*real_part(omega)^6 + 1164.09765419317*imag_part(omega)^24*real_part(omega)^8 + 2793.83437005430*imag_part(omega)^22*real_part(omega)^10 + 5122.02967846394*imag_part(omega)^20*real_part(omega)^12 + 7317.18525490165*imag_part(omega)^18*real_part(omega)^14 + 8231.83341181278*imag_part(omega)^16*real_part(omega)^16 + 7317.18525493145*imag_part(omega)^14*real_part(omega)^18 + 5122.02967846394*imag_part(omega)^12*real_part(omega)^20 + 2793.83437005430*imag_part(omega)^10*real_part(omega)^22 + 1164.09765419317*imag_part(omega)^8*real_part(omega)^24 + 358.183893597743*imag_part(omega)^6*real_part(omega)^26 + ... (huge answer)
Thank you!joaoffFri, 15 Jun 2018 09:23:09 -0500http://ask.sagemath.org/question/42609/convert charpoly()-generated characteristic polynom to symbolic expressionhttp://ask.sagemath.org/question/42486/convert-charpoly-generated-characteristic-polynom-to-symbolic-expression/A characteristic polynom that is generated by the procedure charpoly() has the type:
> sage.rings.polynomial.polynomial_ring.PolynomialRing_field_with_category.element_class.
How can I convert it to one of type 'symbolic expression' without doing 'copy and paste'?
bekalphSun, 03 Jun 2018 12:33:33 -0500http://ask.sagemath.org/question/42486/Polynomial rings with an arbitrary infinite set of variableshttp://ask.sagemath.org/question/42325/polynomial-rings-with-an-arbitrary-infinite-set-of-variables/I have a rational function in $\mathbb{Q}(x)$ which can be computed by substitution of a polynomial in the polynomial ring with integer compositions as variables, with coefficients in $\mathbb{Q}(x)$. The compositions can have unbounded parts. For example, one monomial might be $\frac{1}{1-x}C_{[2,2]}C_{[1,3,1]}$, where $C_\lambda$ is a formal variable related to the integer composition $\lambda$.
I would be happy to have this "symbolic" representation as polynomial over the set of all integer compositions, and not just the polynomial after substituting all of the $C_\lambda$'s. I have tried the following two approaches which failed.
1. Using the `InfinitePolynomialRing` class. The setup used:
<pre><code "python">
sage: R = FractionField(QQ["x"])
sage: R.inject_variables()
Defining x
sage: RC.<c> = InfinitePolynomialRing(R)
sage: f = 1/(1-x)*c[3] + 3*c[2]
</code></pre>
But then f.subs does not work:
sage: f.subs({c[3]:5})
...[snip]
TypeError: <class 'sage.rings.polynomial.infinite_polynomial_ring.InfinitePolynomialGen'> is not hashable
sage: f.subs(c[3]=5) # Expected not to work.
SyntaxError: keyword can't be an expression
There is also a method called `f.specialization` which seems to work with substitution of just one variable, but not with more than one, nor when using an element from `InfinitePolynomialRing(R, "cd")`. There is also the annoyance of using a bijection between the natural numbers and compositions.
2. Using `CombinatorialFreeModule`. I know this might seem as an abuse, but a module can also have the structure of a ring...
<pre><code "python">
sage: R = FractionField(QQ["x"])
sage: FM = FreeAbelianMonoid(Compositions())
sage: CR = Rings().Commutative()
sage: CMR = ModulesWithBasis(R)
sage: M = CombinatorialFreeModule(R, FM, category=(CR,CMR)) # Does not work also with category=None
sage: a = M.an_element(); a*a
...[snip]
TypeError: 'NotImplementedType' object is not callable
</code></pre>
I assume some magic with the `category` argument might help. Note that in my case only a free abelian semigroup indexed by compositions is needed, and not a monoid.
Any help is welcomed.mathzeta2Sun, 13 May 2018 09:10:06 -0500http://ask.sagemath.org/question/42325/cokernel of a map between modules over polynomial ringshttp://ask.sagemath.org/question/42154/cokernel-of-a-map-between-modules-over-polynomial-rings/Define a polynomial ring $R$ as $F_{2}\left[x_{1},x_{1}^{-1},....x_{D},x_{D}^{-1}\right]$ where $D$ is the dimension and $\mathbb{F}_{2}$ is a binary field.
Let $G$ be a free $R$-module of some labels and has rank $t$. $P$
be a free $R$-module of Pauli operators. $\sigma$ is a map from
$G$ to $P$. I want to write a snippet to calculate the cokernel
of this map.
Just for example (taken from page 54 of arxiv.1305.6973 or page 41 of arxiv.1607.01387),
though this is not essential for the question, I can have two ``interaction''
terms in terms of 2-dimensional Pauli operators $X$, $Z$ and Identity
operator $I$ on 4 sites with two 2-dimensional systems per site as
$
II(0,0)-IX(0,1)-XI(1,0)-XX(1,1)
$
and
$
ZZ(0,0)-IZ(0,1)-ZI(1,0)-II(1,1)
$
where on each site $\left(x,y\right)$ (mentioned in the bracket after
the Pauli operators), the first(second) Pauli acts on the first(second)
two dimensional system on that site. The map $\sigma$ can be written
as
$
\sigma=\left(\begin{array}{cc}
y+xy & 0 \\
x+xy & 0 \\
0 & 1+y \\
0 & 1+x
\end{array}\right)
$
where for example, $y+xy$ is a polynomial that specifies the action
on the first two dimensional system as
$
y+x y=0 \hspace{1mm} x^0 y^0+ 1 \hspace{1mm} x^0 y^1+0 \hspace{1mm} x^1 y^0 +1 \hspace{1mm} x^1 y^1
$
where the exponents are the coordinates of the sites and coefficients
$0$ and $1$ imply whether there is a Pauli acting or not.arpitTue, 24 Apr 2018 21:02:38 -0500http://ask.sagemath.org/question/42154/Eigenvalues of matrix with entries in polynomial ringhttp://ask.sagemath.org/question/8797/eigenvalues-of-matrix-with-entries-in-polynomial-ring/Hi!
I just wrote some code on the sage-combinat queue which computes a matrix with entries in a polynomial ring R = PolynomialRing(QQ, 'x', n)
sage: P = Poset(([1,2,3,4], [[1,3],[1,4],[2,3]]), linear_extension = True)
sage: L = P.linear_extensions()
sage: M = L.markov_chain_transition_matrix(labeling = 'source')
sage: M
[-x0 - x1 - x2 x3 x0 + x3 0 0]
[ x1 + x2 -x0 - x1 - x3 0 x1 0]
[ 0 x1 -x0 - x3 0 x1]
[ 0 x0 0 -x0 - x1 - x2 x0 + x3]
[ x0 0 0 x0 + x2 -x0 - x1 - x3]
sage: M.eigenvalues()
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
/Applications/sage-5.0.beta7/devel/sage-combinat/sage/combinat/posets/<ipython console> in <module>()
/Applications/sage-5.0.beta7/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.eigenvalues (sage/matrix/matrix2.c:26415)()
/Applications/sage-5.0.beta7/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.fcp (sage/matrix/matrix2.c:11089)()
/Applications/sage-5.0.beta7/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_element.so in sage.rings.polynomial.polynomial_element.Polynomial.factor (sage/rings/polynomial/polynomial_element.c:22655)()
NotImplementedError:
Is it possible to compute this some other way or is this just not yet implemented (which would surprise me!).
Thanks,
Anne
-----
Edit (originally posted as an answer by the original poster of the question)
Here is an easier example with the question:
sage: R = PolynomialRing(QQ, 'x', 2)
sage: x = R.gens()
sage: M = matrix([[x[0],x[1]],[x[1],x[0]]])
sage: M.eigenvalues()
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
/Applications/sage-5.0.beta7/devel/sage-combinat/sage/combinat/posets/<ipython console> in <module>()
/Applications/sage-5.0.beta7/local/lib/python2.7/site-packages/sage/matrix /matrix2.so in sage.matrix.matrix2.Matrix.eigenvalues (sage/matrix/matrix2.c:26415)()
/Applications/sage-5.0.beta7/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.fcp (sage/matrix/matrix2.c:11089)()
/Applications/sage-5.0.beta7/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_element.so in sage.rings.polynomial.polynomial_element.Polynomial.factor (sage/rings/polynomial/polynomial_element.c:22655)()
NotImplementedError:Anne SchillingFri, 16 Mar 2012 14:05:50 -0500http://ask.sagemath.org/question/8797/Square root of polynomial modulo another irreducible polynomialhttp://ask.sagemath.org/question/42042/square-root-of-polynomial-modulo-another-irreducible-polynomial/Hello,
If I'm not wrong, it is always possible to compute the square root of a polynomial $P$ modulo an irreducible polynomial $g$ when the base field is in $GF(2^m)$, i.e. find $Q \in GF(2^m)$ such that $Q^2 \equiv P \mod g$. Indeed, the operation $Q \rightarrow Q^2 \pmod g$ should be linear (because we are in $GF(2^m)$) so an idea would be to compute the matrix $T$ that perform this operation, and then invert it, but I'd like to find an embedded operation in sage. I tried the sagemath $P.sqrt()$ method, but the problem is that because it does not take into account the modulo, it fails most of the time when the polynomial has some terms with odd power of $X$.
Any idea?
Thanks!tobiasBoraMon, 16 Apr 2018 02:57:45 -0500http://ask.sagemath.org/question/42042/What is the best way to work with ratios of polynomials?http://ask.sagemath.org/question/41760/what-is-the-best-way-to-work-with-ratios-of-polynomials/What is the best way to work with ratios of polynomials with (real and complex) coefficients that can vary from very small to very big values? As I have to keep precision, I am working with symbolic expressions (type 'sage.symbolic.expression.Expression'), but as I operate with the expressions, they become bigger and bigger, and it is more and more difficult to plot them or find its roots, for example.
I perceived that making
# symbolic_expression is some symbolic polynomial with real or complex coefficients
symbolic_expression.polynomial(RR)
symbolic_expression.polynomial(CC)
simplifies the polynomials and eases computations but, as I understand, this converts the symbolic expressions to the real (RR) and complex (CC) rings, rounding the coefficients of the polynomials to machine precision and, in this manner, numerical errors would propagate in my computations. In this scenario, I would need to increase this conversion precision to an arbitrary value that guarantees correct results at the end of all computations.
I've read something about polynomial rings, symbolic ring, etc. I am studying this, but I still don't know what a ring is. Is symbolic ring the same thing as symbolic expressions? Does polynomial rings enable better performance and arbitrary precision?
Cordially
***** EDIT *****
I complement my question to answer to @B r u n o comment.
First, I would like to thank @B r u n o and @slelievre for the tips! I start to better understand things now.
Basically, I generate three polynomials. These polynomials are the basis of my work. I then operate on these polynomials generating other functions, e.g., as the ratio of these polynomials, and continue operating on the new functions, transforming them, plotting, finding the roots of the numerator and denominator polynomials, etc.
I'm slowly improving the code. I started with a plain script, then I created some functions to organize the code, in a procedural approach, and now I am moving to an object-oriented paradigm. This helps me better visualize and understand the problem. I even created a graphical interface to ease the interaction with the program, so as I can easily experiment with different input variables.
The code snippet below reflects the current state of the code it illustrates some inputs to test the code, how I generate my three basic polynomials and some derived functions. Sorry if the code seems silly :-/, I'm not used to develop software, I am still studying.
class FilterPrototype:
def __init(self):
pass
def generate_polynomials(self):
pass
def generate_s21(self, e_poly, p_poly):
return p_poly.polynomial(CC)/e_poly.polynomial(CC)
def generate_s11(self, e_poly, f_poly):
return f_poly.polynomial(CC)/e_poly.polynomial(CC)
def generate_group_delay(self, s21):
return -(1/s21(s=i*omega)*s21(s=i*omega).diff(omega)).imag()
def plot(self, symb_express, bounds):
pl = plot(symb_express, bounds[0], bounds[1])
pl.show()
return list(pl[0])
class ApproxArbitraryPhase(FilterPrototype):
def __init__(self, spec_type, spec, e_order, p_order):
self.spec_type = spec_type # type of specification
self.spec = spec # filter specification
self.e_order = e_order # order of the denominator polynomial
self.p_order = p_order # order of the numerator polynomial
def generate_polynomials(self):
listOmega = [] # List of frequencies (in radians per second)
for a in range(self.e_order + 1):
listOmega.append(a/self.e_order)
if self.spec_type == "group delay":
group_delay_spec(omega) = sage_eval(self.spec, locals={'omega':omega})
elif self.spec_type == "phase":
root.destroy()
sys.exit("You cannot specify a phase. This is yet not implemented.")
else:
root.destroy()
sys.exit("""The spec_type should be "group delay" or "phase".""")
phase(omega) = group_delay_spec(omega).integral(omega)
listPhases = [0]
for a in range(1, self.e_order+1):
listPhases.append(phase(listOmega[a]))
alpha0 = listOmega[1]/tan(listPhases[1])
listAlpha = [alpha0]
for a in range(1, self.e_order):
firstline = listOmega[a+1]^2 - listOmega[a]^2
lastline = alpha0 - listOmega[a+1]/tan(listPhases[a+1])
ic = 1
while ic < a:
lastline = listAlpha[ic] - ((listOmega[a+1]^2 - listOmega[ic]^2)/lastline)
ic+=1
listAlpha.append(firstline/lastline)
ePolynomials = [1, s + listAlpha[0]]
for a in range(2, self.e_order+1):
ePolynomials.append(listAlpha[a-1]*ePolynomials[a-1] + (s^2 + listOmega[a-1]^2)*ePolynomials[a-2])
ePoly = ePolynomials[self.e_order]
pMod = sqrt(ePoly(s=i*omega)*ePoly(s=-i*omega))
integrand = []
for a in range(self.p_order+1):
integrand.append(pMod*chebyshev_T(a, omega)/sqrt(1-omega^2))
c = []
for a in range(self.p_order+1):
if a == 0:
c.append((1/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))
else:
c.append((2/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))
up = 0
for a in range(self.p_order+1):
up = up + c[a]*chebyshev_T(a, omega)
eMod(omega) = (ePoly(s=i*omega)).polynomial(CC)
den = lambda omega: up/abs(eMod(omega))
pl = plot(den(omega), -1, 1) # to improve the precision of norm_factor, the plot_points=200 parameter of the plot function should be augmented
plPoints = []
for pair in list(pl[0]):
plPoints.append(pair[1])
normFactor = Rational(max(plPoints))
pPoly(omega) = up/normFactor
pPoly(s) = pPoly(omega=s/I)
fMod2 = ePoly(s=s)*ePoly(s=-s) - pPoly(s=s)*pPoly(s=-s)
solutions = fMod2.roots(ring=CC)
fPoly = 0
for solution in solutions:
if solution[0].real() < 0 or (solution[0].real() == 0 and solution[0].imag() < 0):
if fPoly == 0:
fPoly = (s - solution[0])
else:
fPoly = fPoly*(s - solution[0])
return (ePoly, fPoly, pPoly)
if __name__ == "__main__":
# execute only if run as a script
s = var("s")
omega = var("omega", domain=RR)
filtro = ApproxArbitraryPhase("group delay", "omega + 5/2", 4, 4)
# filtro = ApproxArbitraryPhase("group delay", "omega + 9/2", 6, 6)
# filtro = ApproxArbitraryPhase("group delay", "-omega + 27/5", 6, 6)
# filtro = ApproxArbitraryPhase("group delay", "piecewise([([0,0.5], 2*omega + 143/20), ((0.5,1), 163/20), ([1,1], 163/20)])", 8, 8)
# filtro = ApproxArbitraryPhase("group delay", "3/5*omega + 106/25", 6, 4)
e,f,p = filtro.generate_polynomials()
print(e)
print(f)
print(p(s))
s21 = filtro.generate_s21(e, p(s))
print(s21)
s11 = filtro.generate_s11(e, f) # possui polinomio do numerador com coeficientes imaginarios
print(s11)
gd = filtro.generate_group_delay(s21)
print(gd)
points_s21 = filtro.plot(20*log(abs(s21(s=i*omega)),10), [-1.5, 1.5])
# print(points_s21)
points_s11 = filtro.plot(20*log(abs(s11(s=i*omega)),10), [-1.5, 1.5])
# print(points_s11)
points_gd = filtro.plot(gd, [-1.5, 1.5])
# print(points_gd)joaoffMon, 26 Mar 2018 04:45:25 -0500http://ask.sagemath.org/question/41760/Dual Quaternion Algebrahttp://ask.sagemath.org/question/41571/dual-quaternion-algebra/Hello,
I'm quite new to SAGE.
I need to work with the algebra of dual quaternions. It can be defined in two equivalent ways:
1. The algebra of quaternions over the dual numbers. Dual numbers are elements of the form a+be where a and b are real and e is such that e^2=0.
2. The algebra of polynomials over the quaternions and variable e modulo e^2.
For the first definition I tried to do this:
P.<e>=PolynomialRing(RR)
S.<e>=P.quo(e*e);
F.<I,J,K> = QuaternionAlgebra(S, -1,-1)
For some reason this doesn't work. What is the best way to define this Dual Quaternion Algebra in SAGE ?
Thank you for your help!TdguerreiroThu, 15 Mar 2018 12:57:36 -0500http://ask.sagemath.org/question/41571/get the coefficients of polynomial of several variables?http://ask.sagemath.org/question/41526/get-the-coefficients-of-polynomial-of-several-variables/ Consider that I have a polynomial with n variables x1,x2,...,xn and want to get the coefficient of that polynomial. For example I can have:
x1+2*x1*x5^2+3*x1^2*x4-5*x2*x5^2+1/2*x2*x3+x6^2*x7*x8+9*x3*x9^3
I would like to ask Sage to give coefficient for the polynomial with variables x1..x9 and the result should be something like this:
coefficient,(exponent of x1..x9);
1,(1,0,0,0,0,0,0,0,0);
2,(1,0,0,0,2,0,0,0,0);
3,(2,0,0,1,0,0,0,0,0);
-5,(0,1,0,0,2,0,0,0,0);
1/2,(0,1,1,0,0,0,0,0,0);
1,(0,0,0,0,0,2,1,1);
9,(0,0,1,0,0,0,0,0,3)
How can I achieve this?
DanialBaghTue, 13 Mar 2018 21:12:52 -0500http://ask.sagemath.org/question/41526/irreducible polynomial defining the finite fieldhttp://ask.sagemath.org/question/41473/irreducible-polynomial-defining-the-finite-field/X = GF(2).polynomial_ring().gen()
F = GF(2^8, name="a", modulus=X^8 + X^6 + X^5 + X + 1)
As example above, the polynomial X^8 + X^6 + X^5 + X + 1 is one of irreducible polynomial defining the finite field GF(2^8).
How can I get all the irreducible polynomial which can define the finite field GF(2^8)?omggggggSat, 10 Mar 2018 11:22:06 -0600http://ask.sagemath.org/question/41473/Factor a quadratic in a quartic polynomialhttp://ask.sagemath.org/question/40304/factor-a-quadratic-in-a-quartic-polynomial/Hi, I have done this calculation using a very tedious way and have checked that it is correct. Can I possibly perform this using Sage only.
I have a polynomial :
D=M^2-(A/(2*p^4))*M+(B/(16*p^4))
where
A=18*p^{10} - 54*p^9 + 59*p^8 + 130*p^7 - 209*p^6 - 98*p^5 + 407*p^4 + 362*p^3 + 49*p^2 - 16*p + 8
and
B=9*( p + 1 )^2*(p^4 - 2*p^3 + 2*p^2 + 2*p + 1)*(4*p^8 - 52*p^7 + 373*p^6 + 68*p^5 - 445*p^4 + 72*p^3 + 163*p^2 - 48*p+ 9).
I have checked using multiple software that the factorization of D using the quartic in $p$
`v^2= p^4-2*p^3+5*p^2+8*p+4` gives
`[M-((A+2*F*v)/(4*p^4))]*[M-((A-2*F*v)/(4*p^4))]` where
`F=9*p^8-18*p^7-7*p^6+45*p^5-21*p^4-74*p^3-18*p^2+6*p-2`.
Can someone help me obtain the same result by using Sage only. Thank you.ShaTue, 26 Dec 2017 09:12:52 -0600http://ask.sagemath.org/question/40304/Solving a system of 18 polynomial equations in sagemathhttp://ask.sagemath.org/question/39998/solving-a-system-of-18-polynomial-equations-in-sagemath/I'm trying to solve a system of 18 polynomial equations in sage.
The system is for sure solvable since I already solved it in Maple and got ALL the possible solutions (only real solutions). In most cases I got between 10 to 25 solutions.
But since Maple is not an open source I cannot really use it.
I'm trying to solve the same equations in sage. but so far I did not succeed.
Since it's the first time I write here, I'm not allowed to attach the Maple and Sage code (two files) that I wanted to attach. Instead, I copied them bellow.
It took Maple 25 minutes to solve the equations (in the code bellow). and it found 22 solutions to the equations (I wrote below all the 22 solutions that Maple found).
I used the command 'Isolate' in Maple to solve the equations (this command is used for solving polynomial equations).
Unfortunately I cannot know how the 'Isolate' command is implemented (I can use the command 'showstat' in Maple and see partial implementation of 'Isolate' since some of the commands that are used inside 'Isolate' are compiled).
But I could see that 'Isolate' uses 'Groebner basis' in order to solve the equations.
Therefore I tried to convert the 18 polynomial equations in sage to groebner basis but the command I.groebner_basis() never finishes (in the example below).
Is there maybe another way of getting the groebner basis of the 18 equations? or any other way of solving the problem (I wrote below) so that the solution is for sure correct?
I'm attaching bellow both Maple's code (that succeeded to solve the equations) and Sage's code (that did not succeed).
I'd appreciate any idea of how to solve the equations.
If it helps, here's the mathematical problem that I'm trying to solve:
given a set of n 3D points P={p_1,p_2,...,p_n} and another set of n 3D lines H={h_1,h_2,...,h_n} (each line h_i is represented by a 3D point and a 3D unit vector). I'm trying to find the rotation matrix R=[r1,r2,r3;r4,r5,r6;r7,r8,r9] and a translation vector T=[t1,t2,t3] such that rotating the set P by R and then translating it by T will result in a new set of points PT={pt_1,...,pt_n} such that the sum of square distances between H and PT (sum_{i=1 to n}(dist(h_i,pt_i)^2)) will be minimal. So this is an optimization problem with constraints (the constraints are R*R^T=I). so there are 12 variables in R and T. Using Lagrange Multipliers algorithm I added 6 more variables b1,...,b6 (the number of constraints from R*R^T=I) and created a new function 'h' with 18 variables. the 18 equations are the derivation of 'h' by all the 18 variables (equal to zero). it is gurantee that one of the solutions to the 18 equations is the (global) minimum for R and T. but for that I have to find ALL the solutions to the equations (22 solutions in this example) just like Maple did.
Thanks
Here is Maple code: (and below sage's code)
restart;
g1:=r1^2+r4^2+r7^2-1;
g2:=r1*r2+r4*r5+r7*r8;
g3:=r1*r3+r4*r6+r7*r9;
g4:=r2^2+r5^2+r8^2-1;
g5:=r2*r3+r5*r6+r8*r9;
g6:=r3^2+r6^2+r9^2-1;
sum_sqr_distances:=(-34.5792590705286*r1+17.0635530183776*r4-2.30671047587914*r7+.533429751140582*t1-18.2777571201152+5.34368522421251*r2-2.6369060119417*r5+.356466130795291*r8-.0824332491501039*t2+31.8951297362284*r3-15.7390369840122*r6+2.12765778880161*r9-.492023587996135*t3)^2+(-63.880273658711*r2+31.5224925463221*r5-4.26132023642126*r8+.98543576110074*t2+5.20415366982094+5.34368522429473*r1-2.63690601164137*r4+.356466130730839*r7-.0824332491771907*t1+5.63520538121055*r3-2.78076015494726*r6+.375912834341308*r9-.0869303242748373*t3)^2+(-31.1892504431195*r3+15.3907123123387*r6-2.08057004842228*r9+.481134487803446*t3+16.460313265992+31.8951297380977*r1-15.7390369788429*r4+2.12765778852522*r7-.492023588007365*t1+5.63520538145412*r2-2.78076015435067*r5+.375912834360444*r8-.086930324248257*t2)^2+(-20.1477543311821*r1+7.55350991666248*r4+20.8920429846501*r7+.506802937949341*t1-13.5125935907312+7.82881106057577*r2-2.9350666574299*r5-8.1180192368411*r8-.19692837123503*t2+18.2686582928059*r3-6.84902591482382*r6-18.9435302828672*r9-.459535566231505*t3)^2+(-36.6286508786965*r2+13.7322935856179*r5+37.9817688885839*r8+.921368583909397*t2-34.3325081149272+7.82881105959312*r1-2.93506665641182*r4-8.11801923407865*r7-.196928371252735*t1+7.29448206284227*r3-2.7347436184352*r6-7.56395131161327*r9-.183487691946913*t3)^2+(-22.7328194787774*r3+8.52266582628087*r6+23.5726043680174*r9+.5718284781815*t3+29.2151845257836+18.2686582944502*r1-6.84902591329222*r4-18.9435302865599*r7-.459535566194367*t1+7.2944820644144*r2-2.73474361877224*r5-7.56395131566163*r8-.183487691915587*t2)^2+(-33.7370053449296*r1+38.9258336170637*r4-60.1132663033064*r7+.996042632734926*t1-4.98803179031169+1.08999592941726*r2-1.25763978603382*r5+1.94217639880877*r8-.0321807583045741*t2+1.8259306385537*r3-2.10676292990297*r6+3.25347948355056*r9-.0539083045951421*t3)^2+(-25.0073507347857*r2+28.8535382444551*r5-44.5585943124607*r8+.738310564652533*t2+4.93804814541239+1.08999592907818*r1-1.25763978578301*r4+1.94217639918338*r7-.0321807583046471*t1+14.8482131236429*r3-17.1319021236029*r6+26.4568411007238*r9-.438374809459876*t3)^2+(-8.9977349393474*r3+10.3816070715198*r6-16.0323428536159*r9+.265646802616713*t3-2.58161819993127+1.82593063830369*r1-2.1067629296129*r4+3.25347948341273*r7-.0539083045952243*t1+14.8482131262288*r2-17.1319021246606*r5+26.4568410945*r8-.438374809459549*t2)^2+(-56.6725654800191*r1-43.1813795970078*r4-58.3887310201578*r7+.980525398272198*t1+49.0916863573367+5.78296846451955*r2+4.40630407967054*r5+5.95808902028457*r8-.100054539766321*t2+5.50887455808305*r3+4.19745958942679*r6+5.67569496876523*r9-.0953122798368383*t3)^2+(-28.0870406673401*r2-21.400781038534*r5-28.937575861748*r8+.48595041543225*t2+19.4018945247715+5.78296846409823*r1+4.40630407900952*r4+5.95808902050728*r7-.100054539765808*t1+28.302910438969*r3+21.5652619383799*r6+29.1599826218892*r9-.489685305322089*t3)^2+(-30.836717169852*r3-23.4958833859599*r6-31.7705184032215*r9+.533524186319051*t3-30.3978529558191+5.50887455772671*r1+4.19745958933706*r4+5.67569496903033*r7-.0953122798304703*t1+28.3029104392002*r2+21.5652619411541*r5+29.1599826221612*r8-.489685305291881*t2)^2+(-57.8205382202822*r1-33.3924384768228*r4-34.894277849845*r7+.850839830317679*t1-12.5672620764049+10.7191476042258*r2+6.19050752309509*r5+6.46892827657879*r8-.157734223996209*t2+21.7070640341081*r3+12.5362340547815*r6+13.1000566054682*r9-.319423430626817*t3)^2+(-56.6217108878986*r2-32.7000933511969*r5-34.1707942023208*r8+.833198866184046*t2-1.62153141633832+10.7191476020505*r1+6.19050752434099*r4+6.46892827744991*r7-.15773422401843*t1+22.9548337765416*r3+13.2568443369017*r6+13.8530766468144*r9-.337784591359295*t3)^2+(-21.4717881919623*r3-12.4003578706288*r6-12.9580693345224*r9+.315961303486945*t3+6.66922254237138+21.7070640313443*r1+12.5362340557455*r4+13.1000566039334*r7-.319423430651098*t1+22.9548337782773*r2+13.2568443352531*r5+13.8530766433259*r8-.337784591337385*t2)^2+(-65.9377723150761*r1-12.678533244879*r4-65.0293030974106*r7+.955841723390729*t1+11.5835324790697+7.68982596884022*r2+1.47860188124137*r5+7.5838780437921*r8-.111472624096221*t2+11.9049207597142*r3+2.28908148248606*r6+11.7408986260292*r9-.17257513524552*t3)^2+(-49.571881779333*r2-9.53169524975009*r5-48.8888964859373*r8+.718599844159693*t2+14.0432100287845+7.6898259693682*r1+1.47860188123288*r4+7.58387804187672*r7-.111472624065637*t1+30.0526392395535*r3+5.77852985094098*r6+29.6385837317152*r9-.435646602436962*t3)^2+(-22.4583184276267*r3-4.31829173040524*r6-22.1488950108346*r9+.32555843242268*t3-12.0350031966363+11.904920759546*r1+2.28908148261245*r4+11.7408986234862*r7-.172575135204147*t1+30.0526392370654*r2+5.77852985129325*r5+29.6385837327811*r8-.435646602452044*t2)^2+(-40.2818933934206*r1+37.3415461911863*r4-17.5121421408269*r7+.688992135032061*t1+41.1079956272437+18.3480006638221*r2-17.0087018409058*r5+7.97660607993154*r8-.313829045441664*t2+19.8946971640516*r3-18.4424983689296*r6+8.64901659965967*r9-.340284150555745*t3)^2+(-39.9505168284587*r2+37.0343582140755*r5-17.3680795584234*r8+.68332418286889*t2-19.0229410779165+18.3480006624515*r1-17.0087018386759*r4+7.97660607742464*r7-.313829045473842*t1+20.0751637635565*r3-18.6097919416627*r6+8.72747261243189*r9-.343370898899277*t3)^2+(-36.6974975173012*r3+34.0187906543919*r6-15.9538626068613*r9+.627683682093741*t3-20.0272581058992+19.8946971612448*r1-18.4424983632605*r4+8.64901659883117*r7-.340284150571336*t1+20.0751637622239*r2-18.609791938382*r5+8.72747261433876*r8-.343370898879803*t2)^2+(-19.9003018252227*r1+2.10723968998899*r4-5.99083386660634*r7+.268355364500816*t1-6.36217996635187+32.4452671740533*r2-3.43562400998522*r5+9.76739986562233*r8-.437524092703596*t2+5.19806634264006*r3-.550422390906632*r6+1.56483817002714*r9-.0700958709468267*t3)^2+(-54.7542305404035*r2+5.79791647536792*r5-16.4833428943737*r8+.738360233327084*t2+12.3293109763825+32.4452671726755*r1-3.43562401018043*r4+9.76739986673756*r7-.437524092757061*t1+3.1084479408544*r3-.329153041695288*r6+.935774510511675*r9-.0419173883795401*t3)^2+(-73.6585215390464*r3+7.79968874305612*r6-22.1743353083965*r9+.99328440217389*t3-10.5500615738974+5.19806634229016*r1-.550422390946524*r4+1.5648381697327*r7-.0700958709446797*t1+3.10844794077716*r2-.329153041700441*r5+.935774510228754*r8-.041917388373134*t2)^2;
h:=sum_sqr_distances+g1*b1+g2*b2+g3*b3+g4*b4+g5*b5+g6*b6;
diff_t1:=diff(h,t1);
diff_t2:=diff(h,t2);
diff_t3:=diff(h,t3);
diff_r1:=diff(h,r1);
diff_r2:=diff(h,r2);
diff_r3:=diff(h,r3);
diff_r4:=diff(h,r4);
diff_r5:=diff(h,r5);
diff_r6:=diff(h,r6);
diff_r7:=diff(h,r7);
diff_r8:=diff(h,r8);
diff_r9:=diff(h,r9);
diff_b1:=diff(h,b1);
diff_b2:=diff(h,b2);
diff_b3:=diff(h,b3);
diff_b4:=diff(h,b4);
diff_b5:=diff(h,b5);
diff_b6:=diff(h,b6);
vars := [op(indets(h, And(name, Non(constant))))];
polysys:={diff_t1,diff_t2,diff_t3,diff_r1,diff_r4,diff_r7,diff_r2,diff_r5,diff_r8,diff_r3,diff_r6,diff_r9,diff_b1,diff_b2,diff_b3,diff_b4,diff_b5,diff_b6};
sols := CodeTools:-Usage(RootFinding:-Isolate(polysys, vars, output = numeric, method = RS));
And this is sage code (that couldn't find the groebner basis):
P.<r1,r2,r3,r4,r5,r6,r7,r8,r9,t1,t2,t3,b1,b2,b3,b4,b5,b6>=PolynomialRing(QQ,order='degrevlex')
g1=r1^2+r4^2+r7^2-1
g2=r1*r2+r4*r5+r7*r8
g3=r1*r3+r4*r6+r7*r9
g4=r2^2+r5^2+r8^2-1
g5=r2*r3+r5*r6+r8*r9
g6=r3^2+r6^2+r9^2-1
sum_sqr_distances=(-34.5792590705286*r1+17.0635530183776*r4-2.30671047587914*r7+.533429751140582*t1-18.2777571201152+5.34368522421251*r2-2.6369060119417*r5+.356466130795291*r8-.0824332491501039*t2+31.8951297362284*r3-15.7390369840122*r6+2.12765778880161*r9-.492023587996135*t3)^2+(-63.880273658711*r2+31.5224925463221*r5-4.26132023642126*r8+.98543576110074*t2+5.20415366982094+5.34368522429473*r1-2.63690601164137*r4+.356466130730839*r7-.0824332491771907*t1+5.63520538121055*r3-2.78076015494726*r6+.375912834341308*r9-.0869303242748373*t3)^2+(-31.1892504431195*r3+15.3907123123387*r6-2.08057004842228*r9+.481134487803446*t3+16.460313265992+31.8951297380977*r1-15.7390369788429*r4+2.12765778852522*r7-.492023588007365*t1+5.63520538145412*r2-2.78076015435067*r5+.375912834360444*r8-.086930324248257*t2)^2+(-20.1477543311821*r1+7.55350991666248*r4+20.8920429846501*r7+.506802937949341*t1-13.5125935907312+7.82881106057577*r2-2.9350666574299*r5-8.1180192368411*r8-.19692837123503*t2+18.2686582928059*r3-6.84902591482382*r6-18.9435302828672*r9-.459535566231505*t3)^2+(-36.6286508786965*r2+13.7322935856179*r5+37.9817688885839*r8+.921368583909397*t2-34.3325081149272+7.82881105959312*r1-2.93506665641182*r4-8.11801923407865*r7-.196928371252735*t1+7.29448206284227*r3-2.7347436184352*r6-7.56395131161327*r9-.183487691946913*t3)^2+(-22.7328194787774*r3+8.52266582628087*r6+23.5726043680174*r9+.5718284781815*t3+29.2151845257836+18.2686582944502*r1-6.84902591329222*r4-18.9435302865599*r7-.459535566194367*t1+7.2944820644144*r2-2.73474361877224*r5-7.56395131566163*r8-.183487691915587*t2)^2+(-33.7370053449296*r1+38.9258336170637*r4-60.1132663033064*r7+.996042632734926*t1-4.98803179031169+1.08999592941726*r2-1.25763978603382*r5+1.94217639880877*r8-.0321807583045741*t2+1.8259306385537*r3-2.10676292990297*r6+3.25347948355056*r9-.0539083045951421*t3)^2+(-25.0073507347857*r2+28.8535382444551*r5-44.5585943124607*r8+.738310564652533*t2+4.93804814541239+1.08999592907818*r1-1.25763978578301*r4+1.94217639918338*r7-.0321807583046471*t1+14.8482131236429*r3-17.1319021236029*r6+26.4568411007238*r9-.438374809459876*t3)^2+(-8.9977349393474*r3+10.3816070715198*r6-16.0323428536159*r9+.265646802616713*t3-2.58161819993127+1.82593063830369*r1-2.1067629296129*r4+3.25347948341273*r7-.0539083045952243*t1+14.8482131262288*r2-17.1319021246606*r5+26.4568410945*r8-.438374809459549*t2)^2+(-56.6725654800191*r1-43.1813795970078*r4-58.3887310201578*r7+.980525398272198*t1+49.0916863573367+5.78296846451955*r2+4.40630407967054*r5+5.95808902028457*r8-.100054539766321*t2+5.50887455808305*r3+4.19745958942679*r6+5.67569496876523*r9-.0953122798368383*t3)^2+(-28.0870406673401*r2-21.400781038534*r5-28.937575861748*r8+.48595041543225*t2+19.4018945247715+5.78296846409823*r1+4.40630407900952*r4+5.95808902050728*r7-.100054539765808*t1+28.302910438969*r3+21.5652619383799*r6+29.1599826218892*r9-.489685305322089*t3)^2+(-30.836717169852*r3-23.4958833859599*r6-31.7705184032215*r9+.533524186319051*t3-30.3978529558191+5.50887455772671*r1+4.19745958933706*r4+5.67569496903033*r7-.0953122798304703*t1+28.3029104392002*r2+21.5652619411541*r5+29.1599826221612*r8-.489685305291881*t2)^2+(-57.8205382202822*r1-33.3924384768228*r4-34.894277849845*r7+.850839830317679*t1-12.5672620764049+10.7191476042258*r2+6.19050752309509*r5+6.46892827657879*r8-.157734223996209*t2+21.7070640341081*r3+12.5362340547815*r6+13.1000566054682*r9-.319423430626817*t3)^2+(-56.6217108878986*r2-32.7000933511969*r5-34.1707942023208*r8+.833198866184046*t2-1.62153141633832+10.7191476020505*r1+6.19050752434099*r4+6.46892827744991*r7-.15773422401843*t1+22.9548337765416*r3+13.2568443369017*r6+13.8530766468144*r9-.337784591359295*t3)^2+(-21.4717881919623*r3-12.4003578706288*r6-12.9580693345224*r9+.315961303486945*t3+6.66922254237138+21.7070640313443*r1+12.5362340557455*r4+13.1000566039334*r7-.319423430651098*t1+22.9548337782773*r2+13.2568443352531*r5+13.8530766433259*r8-.337784591337385*t2)^2+(-65.9377723150761*r1-12.678533244879*r4-65.0293030974106*r7+.955841723390729*t1+11.5835324790697+7.68982596884022*r2+1.47860188124137*r5+7.5838780437921*r8-.111472624096221*t2+11.9049207597142*r3+2.28908148248606*r6+11.7408986260292*r9-.17257513524552*t3)^2+(-49.571881779333*r2-9.53169524975009*r5-48.8888964859373*r8+.718599844159693*t2+14.0432100287845+7.6898259693682*r1+1.47860188123288*r4+7.58387804187672*r7-.111472624065637*t1+30.0526392395535*r3+5.77852985094098*r6+29.6385837317152*r9-.435646602436962*t3)^2+(-22.4583184276267*r3-4.31829173040524*r6-22.1488950108346*r9+.32555843242268*t3-12.0350031966363+11.904920759546*r1+2.28908148261245*r4+11.7408986234862*r7-.172575135204147*t1+30.0526392370654*r2+5.77852985129325*r5+29.6385837327811*r8-.435646602452044*t2)^2+(-40.2818933934206*r1+37.3415461911863*r4-17.5121421408269*r7+.688992135032061*t1+41.1079956272437+18.3480006638221*r2-17.0087018409058*r5+7.97660607993154*r8-.313829045441664*t2+19.8946971640516*r3-18.4424983689296*r6+8.64901659965967*r9-.340284150555745*t3)^2+(-39.9505168284587*r2+37.0343582140755*r5-17.3680795584234*r8+.68332418286889*t2-19.0229410779165+18.3480006624515*r1-17.0087018386759*r4+7.97660607742464*r7-.313829045473842*t1+20.0751637635565*r3-18.6097919416627*r6+8.72747261243189*r9-.343370898899277*t3)^2+(-36.6974975173012*r3+34.0187906543919*r6-15.9538626068613*r9+.627683682093741*t3-20.0272581058992+19.8946971612448*r1-18.4424983632605*r4+8.64901659883117*r7-.340284150571336*t1+20.0751637622239*r2-18.609791938382*r5+8.72747261433876*r8-.343370898879803*t2)^2+(-19.9003018252227*r1+2.10723968998899*r4-5.99083386660634*r7+.268355364500816*t1-6.36217996635187+32.4452671740533*r2-3.43562400998522*r5+9.76739986562233*r8-.437524092703596*t2+5.19806634264006*r3-.550422390906632*r6+1.56483817002714*r9-.0700958709468267*t3)^2+(-54.7542305404035*r2+5.79791647536792*r5-16.4833428943737*r8+.738360233327084*t2+12.3293109763825+32.4452671726755*r1-3.43562401018043*r4+9.76739986673756*r7-.437524092757061*t1+3.1084479408544*r3-.329153041695288*r6+.935774510511675*r9-.0419173883795401*t3)^2+(-73.6585215390464*r3+7.79968874305612*r6-22.1743353083965*r9+.99328440217389*t3-10.5500615738974+5.19806634229016*r1-.550422390946524*r4+1.5648381697327*r7-.0700958709446797*t1+3.10844794077716*r2-.329153041700441*r5+.935774510228754*r8-.041917388373134*t2)^2
sum_sqr_distances=P(sum_sqr_distances)
h=sum_sqr_distances+g1*b1+g2*b2+g3*b3+g4*b4+g5*b5+g6*b6
diff_t1=diff(h,t1)
diff_t2=diff(h,t2)
diff_t3=diff(h,t3)
diff_r1=diff(h,r1)
diff_r2=diff(h,r2)
diff_r3=diff(h,r3)
diff_r4=diff(h,r4)
diff_r5=diff(h,r5)
diff_r6=diff(h,r6)
diff_r7=diff(h,r7)
diff_r8=diff(h,r8)
diff_r9=diff(h,r9)
diff_b1=diff(h,b1)
diff_b2=diff(h,b2)
diff_b3=diff(h,b3)
diff_b4=diff(h,b4)
diff_b5=diff(h,b5)
diff_b6=diff(h,b6)
I=Ideal(diff_t1,diff_t2,diff_t3,diff_r1,diff_r2,diff_r3,diff_r4,diff_r5,diff_r6,diff_r7,diff_r8,diff_r9,diff_b1,diff_b2,diff_b3,diff_b4,diff_b5,diff_b6)
I.groebner_basis()
I.variety(RR)
And these are the 22 solutions the Maple found (all the real solutions that solve the equations):
The goal is to get the same solutions from sage.
[[b1 = -6405.651236, b2 = 1097.011411, b3 = 6432.603484, b4 = -3486.767199, b5 = 2786.505511, b6 = -2066.461540, r1 = -.2735272365, r2 = -.4584555697, r3 = -.8455775195, r4 = .9594115168, r5 = -0.6729917007e-1, r6 = -.2738619418, r7 = -0.6864686727e-1, r8 = .8861655107, r9 = -.4582557095, t1 = -31.53181188, t2 = -3.035192389, t3 = -64.33835896],
[b1 = -663.0587184, b2 = 1206.645590, b3 = -935.8274112, b4 = -2827.613858, b5 = 1260.996624, b6 = -676.9562322, r1 = .8750945364, r2 = -.1796187040, r3 = -.4493847723, r4 = -.4474525861, r5 = -.6540696330, r6 = -.6099008922, r7 = -.1843793253, r8 = .7347993170, r9 = -.6527436159, t1 = 24.66279957, t2 = 6.035593829, t3 = -55.14919482],
[b1 = -5334.007938, b2 = -147.4788994, b3 = 6381.507728, b4 = -1654.278431, b5 = 3038.269718, b6 = -2677.888758, r1 = 0.1140465611e-1, r2 = -.8980265159, r3 = -.4397934863, r4 = .8207951804, r5 = .2596099997, r6 = -.5088201253, r7 = -.5711087512, r8 = .3551774554, r9 = -.7400565989, t1 = -38.50311082, t2 = -47.94244185, t3 = -52.00642291],
[b1 = -3711.964527, b2 = 3878.643381, b3 = -871.0350173, b4 = -4720.758037, b5 = 5739.216423, b6 = -3328.445679, r1 = -.3452025304, r2 = -.5900255619, r3 = -.7298664599, r4 = -.5720856277, r5 = .7487785844, r6 = -.3347367116, r7 = .7440115910, r8 = .3019941520, r9 = -.5960254060, t1 = 4.420063381, t2 = -25.83338983, t3 = -49.39760950],
[b1 = -6771.196603, b2 = 3126.930183, b3 = 3836.296676, b4 = -2119.174006, b5 = -1053.266216, b6 = -647.0603393, r1 = -0.5288779601e-1, r2 = .2948075346, r3 = -.9540919235, r4 = .9579998211, r5 = -.2546859938, r6 = -.1318005592, r7 = .2818496477, r8 = .9209905331, r9 = .2689557845, t1 = -1.568329794, t2 = 41.81946440, t3 = -47.91552928],
[b1 = -426.4877331, b2 = 1777.234977, b3 = -1396.594665, b4 = -4704.976941, b5 = 6548.202362, b6 = -2797.527674, r1 = .6772160397, r2 = -.2551001296, r3 = -.6901466217, r4 = -.3808728417, r5 = .6809888166, r6 = -.6254519247, r7 = .6295349982, r8 = .6864241804, r9 = .3640158382, t1 = 51.90734298, t2 = -10.61796786, t3 = -36.65155778],
[b1 = -8345.957097, b2 = -5992.141964, b3 = 9771.580880, b4 = -9281.942970, b5 = 14067.10170, b6 = -7238.069438, r1 = -.4112595863, r2 = -.5902847582, r3 = -.6945714196, r4 = -.9059223617, r5 = .1803855392, r6 = .3831001589, r7 = .1008475446, r8 = -.7867813937, r9 = .6089374444, t1 = -25.97531493, t2 = -66.23689116, t3 = -31.46159720],
[b1 = -5343.147703, b2 = 183.2763579, b3 = 5426.240402, b4 = -11788.39198, b5 = 16609.81042, b6 = -7181.846374, r1 = -.2988639430, r2 = -.6827291572, r3 = -.6667542587, r4 = .7307518828, r5 = -.6130817579, r6 = .3002206588, r7 = -.6137442703, r8 = -.3975068000, r9 = .6821336486, t1 = -51.64655336, t2 = -55.45375783, t3 = -29.08482105],
[b1 = -4075.023103, b2 = 6883.346857, b3 = -1945.068213, b4 = -3311.752332, b5 = 3039.918275, b6 = -921.4062369, r1 = .5825663050, r2 = .8122102291, r3 = 0.3051301418e-1, r4 = -.8079724338, r5 = .5827884450, r6 = -0.8682266235e-1, r7 = 0.8830088656e-1, r8 = -0.2592628327e-1, r9 = -.9957563865, t1 = 23.77615250, t2 = 40.91304243, t3 = -28.03945694],
[b1 = -6187.757714, b2 = 838.7708899, b3 = 6066.514065, b4 = -441.1250637, b5 = -851.5589307, b6 = -1128.553632, r1 = .2065219306, r2 = .9631666887, r3 = -.1722167877, r4 = .8113374284, r5 = -.2669555392, r6 = -.5200637629, r7 = -.5468823178, r8 = -0.3232135326e-1, r9 = -.8365853576, t1 = -21.72396020, t2 = 62.90274148, t3 = -22.56105906],
[b1 = -3971.638517, b2 = 5211.591913, b3 = -2909.389068, b4 = -3415.632999, b5 = 2348.040874, b6 = -1624.438321, r1 = -.6376325216, r2 = -.7577247398, r3 = -.1388451874, r4 = -.6573617358, r5 = .6291742407, r6 = -.4147473004, r7 = .4016221056, r8 = -.1731848536, r9 = -.8992812078, t1 = -24.08083201, t2 = -42.77528618, t3 = -16.88667732],
[b1 = -1107.397052, b2 = -1352.354473, b3 = 559.8732630, b4 = -682.7595891, b5 = 1919.821403, b6 = -459.8284277, r1 = -.5540008534, r2 = .8147846552, r3 = -.1709064659, r4 = .3838720628, r5 = .4321713781, r6 = .8160086638, r7 = .7387322206, r8 = .3864632785, r9 = -.5521963786, t1 = -3.462811975, t2 = 66.31134317, t3 = -11.39008932],
[b1 = -1858.210974, b2 = -1443.942073, b3 = -117.7029385, b4 = -541.5565807, b5 = 635.7869277, b6 = -1292.459322, r1 = -.6056473654, r2 = .7430910328, r3 = -.2846172619, r4 = -.4570833521, r5 = -.6176622647, r6 = -.6399751058, r7 = .6513571049, r8 = .2575054246, r9 = -.7137400635, t1 = -8.170594439, t2 = 73.18220544, t3 = -5.574473809],
[b1 = -2879.955409, b2 = -7247.457464, b3 = 7603.360417, b4 = -10386.60081, b5 = 15680.46439, b6 = -7326.146084, r1 = .8898800024, r2 = -.3777475308, r3 = .2557740885, r4 = -.2159488982, r5 = .1450669474, r6 = .9655680474, r7 = -.4018453119, r8 = -.9144738289, r9 = 0.4751801242e-1, t1 = 25.17284844, t2 = -56.81452504, t3 = 4.483386705],
[b1 = -2594.500775, b2 = -8451.746216, b3 = 8378.142403, b4 = -8445.314669, b5 = 16539.21636, b6 = -7726.257147, r1 = -.8887915121, r2 = .4283658935, r3 = .1629487933, r4 = .2014992422, r5 = 0.4589445400e-1, r6 = .9784128753, r7 = -.4116402597, r8 = -.9024391172, r9 = .1271060041, t1 = -72.41708628, t2 = -1.180372243, t3 = 11.62696253],
[b1 = -3386.721976, b2 = -9232.755155, b3 = 8735.844228, b4 = -8449.692294, b5 = 14675.00399, b6 = -7470.241638, r1 = .7020126763, r2 = -.7115667558, r3 = -0.2917112402e-1, r4 = -.6916841933, r5 = -.6714995329, r6 = -.2658220344, r7 = -.1695617265, r8 = -.2067876432, r9 = .9635806617, t1 = 22.56779868, t2 = -57.27997094, t3 = 12.66120812],
[b1 = -2344.035915, b2 = 3947.708934, b3 = -879.5603548, b4 = -4851.200663, b5 = 7270.443175, b6 = -2742.299953, r1 = .4417879292, r2 = .7119951524, r3 = .5457896376, r4 = -.5754568251, r5 = .6916351055, r6 = -.4364519713, r7 = .6882389614, r8 = .1212591594, r9 = -.7152785110, t1 = 49.50818048, t2 = 50.77308450, t3 = 32.05737927],
[b1 = -4793.592445, b2 = -9508.057512, b3 = 9867.840732, b4 = -7363.405470, b5 = 13671.66420, b6 = -8641.244153, r1 = -.8607833310, r2 = .4979577238, r3 = .1053098401, r4 = -.5071531704, r5 = -.8216734210, r6 = -.2600931582, r7 = -0.4298510042e-1, r8 = -.2772920743, r9 = .9598236227, t1 = -55.37291300, t2 = 22.92886207, t3 = 35.86152727],
[b1 = -1143.596369, b2 = 1527.374908, b3 = 1286.733075, b4 = -1246.223196, b5 = 880.0474364, b6 = -176.7812153, r1 = .5096005566, r2 = -.5301666308, r3 = .6776655638, r4 = -0.1872847237e-1, r5 = .7805852008, r6 = .6247687481, r7 = .8602072524, r8 = .3310741426, r9 = -.3878574417, t1 = 62.72704338, t2 = -19.19220727, t3 = 40.08281887],
[b1 = -6277.680746, b2 = -547.8928701, b3 = 6669.122409, b4 = -11168.08424, b5 = 15793.41027, b6 = -8033.605502, r1 = .4599469977, r2 = .3665077232, r3 = .8087773786, r4 = .5464600340, r5 = -.8347594860, r6 = 0.6751319714e-1, r7 = -.6998786970, r8 = -.4109120215, r9 = .5842269422, t1 = -10.38583601, t2 = 10.41107518, t3 = 61.92902662],
[b1 = -9151.429703, b2 = -6008.117168, b3 = 11049.83952, b4 = -7190.606039, b5 = 13545.27447, b6 = -7223.768017, r1 = 0.7425420274e-1, r2 = .5809660735, r3 = .8105336112, r4 = -.9944330906, r5 = .1040687725, r6 = 0.1650814969e-1, r7 = -0.7476056307e-1, r8 = -.8072472434, r9 = .5854594317, t1 = -4.663325049, t2 = 4.554915592, t3 = 62.58462359],
[b1 = -1226.728760, b2 = 1147.011945, b3 = -400.7753802, b4 = -4507.321778, b5 = 8005.881456, b6 = -3241.483736, r1 = -0.9032875350e-2, r2 = .7488743127, r3 = .6626504893, r4 = -.2759411564, r5 = .6350798638, r6 = -.7214776814, r7 = .9611320853, r8 = .1893695603, r9 = -.2009086466, t1 = 40.46849908, t2 = 60.56033489, t3 = 65.08406719]]david_cThu, 07 Dec 2017 12:00:21 -0600http://ask.sagemath.org/question/39998/Change Precision on complex_roots()http://ask.sagemath.org/question/39988/change-precision-on-complex_roots/ I am trying to find the complex roots of the polynomial
poly = x^7 - 6*x^6 + 15*x^5 - 20*x^4 + 15*x^3 - 6*x^2 + x
But when I do poly.complex_roots(), the system gives:
PariError: overflow in expo()
Apparently there are options for how much precision you want when computing roots -- one option is to use Pari, which is the high-precision option, and the other NumPy, which is the low-precision option. The default is set to use Pari, which apparently overloads when I try to compute the roots of this polynomial (and many others as well, this polynomial is just one example).
How do I change the complex_roots() function to get lower-precision roots?
Alternatively, how do I deal with the PariError?
cshiringWed, 06 Dec 2017 23:13:28 -0600http://ask.sagemath.org/question/39988/How do I code a Laurent Series with variable coefficients?http://ask.sagemath.org/question/39662/how-do-i-code-a-laurent-series-with-variable-coefficients/Edit: Adding more context.
I am attempting the following procedure:
1. Begin with a polynomial $Z(u)$ with variable coefficients, of the form $1 + a*u + b*u^2 + c*u^3 + p*b*u^4 + p^2*a*u^5 + p^3*u^6$.
2. Examine the coefficients of $Z'(u)/Z(u)$ as a power series in $u$.
It is this quest which leads me to attempt to construct a LaurentSeriesRing with variable coefficients. However, I keep encountering TypeErrors, I am wondering if a kind soul could help me in my quest. I will use here a very simple polynomial f to get the point across.
I am attempting to construct a LaurentSeriesRing with variable coefficients. However, I keep encountering TypeErrors, I am wondering if a kind soul could help me in my quest. I will use here a very simple polynomial f to get the point across.
sage: R.<t> = LaurentSeriesRing(QQ, 't')
sage: var('a')
a
sage: f = 1 + a*t
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-31-06d3f2f41e45> in <module>()
----> 1 f = Integer(1) + a*t
sage/structure/element.pyx in sage.structure.element.Element.__mul__ (/usr/lib/sagemath//src/build/cythonized/sage/structure/element.c:12443)()
sage/structure/coerce.pyx in sage.structure.coerce.CoercionModel_cache_maps.bin_op (/usr/lib/sagemath//src/build/cythonized/sage/structure/coerce.c:10496)()
TypeError: unsupported operand parent(s) for '*': 'Symbolic Ring' and 'Laurent Series Ring in t over Rational Field'
I also tried:
sage: R.<u> = QQ[]
sage: var('a')
a
sage: f = 1 + a*u
sage: ff = derivative(f, u)
sage: R.<u> = LaurentSeriesRing(QQ); R
sage: f/ff + O(u^5)
----------------------
TypeError Traceback (most recent call last)
<ipython-input-28-c4846de7ced8> in <module>()
----> 1 f/ff + O(u**Integer(5))
sage/structure/element.pyx in sage.structure.element.Element.__add__ (/usr/lib/sagemath//src/build/cythonized/sage/structure/element.c:11198)()
sage/structure/coerce.pyx in sage.structure.coerce.CoercionModel_cache_maps.bin_op (/usr/lib/sagemath//src/build/cythonized/sage/structure/coerce.c:10496)()
TypeError: unsupported operand parent(s) for '+': 'Symbolic Ring' and 'Laurent Series Ring in u over Rational Field'masseygirlSun, 19 Nov 2017 13:35:21 -0600http://ask.sagemath.org/question/39662/Factoring polynomials with symbolic expressionshttp://ask.sagemath.org/question/39772/factoring-polynomials-with-symbolic-expressions/ This question is about writing some code that factors a desired polynomial within $\mathbb{Q}$, $\mathbb{R}$ and $\mathbb{C}$ for educational purposes.
So far, I have the following code:
@interact
def _(p = input_box(default=x^5 + x^4 - 8*x^3 + 11*x^2 - 15*x + 2, label = 'P(x) = ')):
q.<x> = PolynomialRing(QQ, 'x')
r.<x> = PolynomialRing(RR, 'x')
c.<x> = PolynomialRing(CC, 'x')
a.<x> = PolynomialRing(AA, 'x')
qp = PolynomialRing(QQ, 'x')(p)
print(factor(qp))
rp = r(p)
print(factor(rp))
cp = c(p)
print(factor(cp))
ap = a(p)
print(factor(ap))
whose output is
(x - 2) * (x^4 + 3*x^3 - 2*x^2 + 7*x - 1)
(x - 2.00000000000000) * (x - 0.147637797932293) * (x + 3.96552349222940) * (x^2 - 0.817885694297105*x + 1.70805524786870)
(x - 2.00000000000000) * (x - 0.408942847148552 - 1.24129810909174*I) * (x - 0.408942847148552 + 1.24129810909174*I) * (x - 0.147637797932293) * (x + 3.96552349222940)
(x - 2.000000000000000?) * (x - 0.1476377979322930?) * (x + 3.965523492229398?) * (x^2 - 0.8178856942971048?*x + 1.708055247868697?)
The factorization is perfect, but I would like the output to be symbolic whenever possible (ie. written with radicals instead of decimals). I know about the use of the ring 'AA' and that the numbers ending with '?' may be translated into symbolic, but I do not know the smartest way to do this.
Which would be the smartest way to achieve this goal? In the same way, ideas to improve the previous code are also welcomejepstraFri, 24 Nov 2017 11:17:56 -0600http://ask.sagemath.org/question/39772/How can I readline a one-term polynomial from a user?http://ask.sagemath.org/question/39760/how-can-i-readline-a-one-term-polynomial-from-a-user/I would like to ask the user for a function. How can I do that? Is there something like readline?veritasFri, 24 Nov 2017 02:08:11 -0600http://ask.sagemath.org/question/39760/Characteristic polynomial wont be used in solvehttp://ask.sagemath.org/question/39746/characteristic-polynomial-wont-be-used-in-solve/ When I try to find roots in a characteristic polynomial it gives me errors:
sage: #Diagonalmatrix
....:
....:
....: A=matrix([[1,-1,2],
....: [-1,1,2],
....: [2,2,-2]])
....: var('x')
....: poly=A.characteristic_polynomial()
....: eq1=solve(poly==0,x)
....:
x
--------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-86-fee7a7de2ea1> in <module>()
7 var('x')
8 poly=A.characteristic_polynomial()
----> 9 eq1=solve(poly==Integer(0),x)
/usr/lib/python2.7/site-packages/sage/symbolic/relation.pyc in solve(f, *args, **kwds)
816
817 if not isinstance(f, (list, tuple)):
--> 818 raise TypeError("The first argument must be a symbolic expression or a list of symbolic expressions.")
819
820 if len(f)==1:
TypeError: The first argument must be a symbolic expression or a list of symbolic expressions.
sage:
What can I do in order to use the polynomial in an equation I wish to solve?PoetastropheThu, 23 Nov 2017 11:18:11 -0600http://ask.sagemath.org/question/39746/Random matrix satisfying a given polynomialhttp://ask.sagemath.org/question/39721/random-matrix-satisfying-a-given-polynomial/ If a polynomial f(x) of order n is given, can we find a random square matrix A of order m so that f(A)=0?
I tried to construct it by finding the roots of f(x) and then creating random matrix with those roots as eigenvalues. But the problem occurs when n is not equal to m. I'm unable to set the eigenvalue, dimensions suitably.Deepak SarmaWed, 22 Nov 2017 04:38:13 -0600http://ask.sagemath.org/question/39721/Segmentation fault when multiplying by variablehttp://ask.sagemath.org/question/39204/segmentation-fault-when-multiplying-by-variable/In Sage 7.5.1 I'm trying to work with unknown values in GF(3) and polynomials.
The following code gives me a *segmentation fault*:
P.<x> = GF(3)['x']
var('a')
sigma = 2*x+ 1
print(a*(sigma))
What is the proper way to handle unknown GF(3) values like "a" in Sage?PsiWed, 18 Oct 2017 15:30:48 -0500http://ask.sagemath.org/question/39204/Addition polynomial in finite field errorhttp://ask.sagemath.org/question/39195/addition-polynomial-in-finite-field-error/I try to do some additions with finite field (GF) just like below. What happen with A?
R.<x> = GF(2)[]
K.<x> = GF(2^8, modulus=x^8 + x^4 + x^3 + x + 1)
A = 1 + 1
B = x + x
C = x + x^2
print A
print B
print C
Result:
2
0
x^2 + xNiyamabrataTue, 17 Oct 2017 09:46:19 -0500http://ask.sagemath.org/question/39195/Is there a way to get the homogeneous part of certain degree of a (multivariate) polynomial?http://ask.sagemath.org/question/39177/is-there-a-way-to-get-the-homogeneous-part-of-certain-degree-of-a-multivariate-polynomial/Every multivariate polynomial $f\in\Bbb F[x_1,\ldots,x_n]$ of degree $d$ can be written as $f = f_0+f_1+\cdots+f_d$, where $f_i$ is a **homogeneous** polynomial of degree $i$. Is there a *direct* way to get each $f_i$ given $f$ in SageMath? For a specific application I have where I only need $f_d$ I am homogenizing and then setting $h=0$, and based on this I wrote an ugly script that recursively finds $f_i$.
Is there a cleaner (and more efficient) way to do this?
Thanks for the help!
**EDIT:**
This is the code I'm using to obtain $f_d$ from $f$
fd = R( f.homogenize()(h=0) )
where R is the multivariate polynomial ring (parent of $f$). If I want $f_{d-1}$ for example, I can define $g$ as $f - f_d$ and apply the line to $g$. This recursive definition is not satisfactory since to get $f_i$ I need to have all $f_{i+1},\ldots,f_d$ first, which is inefficient. Also, that trick of homogenizing, evaluating $h=0$ and coercing the result back to the original polynomial ring is not very neat.descuderoSun, 15 Oct 2017 12:37:26 -0500http://ask.sagemath.org/question/39177/