expressing f(x)/g(x) as a polynomial

68 views
Skip to first unread message

gosia

unread,
Dec 22, 2017, 5:27:06 PM12/22/17
to sympy
Dear SymPy users,

I'm new to SymPy and i'm struggling for some time now with expressing a polynomial fraction as a polynomial.

I try to calculate:
w(x) = f(x^10)/g(x^2)
which should be a function of x (of the order 8)

what i get (using various sympy functions) is:
w(x) = a_8 * x^8 + ... + a_0 + p(x)/q(x^2)
or
w(x) = a_8 * x^8 + ... + a_0 + p(x^0)/q(x^2)

instead of:
w(x) = a_8 * x^8 + ... + a_0


below is the full story with a part of a script i'm using (i installed sympy 1.1.1 with python 3.6 in anaconda environment; i'm using linux ubuntu 16.04)

any hints will be extremely helpful!
best,
gosia


------------------------------------------------------------------------
import sympy as sp

# 1. i start with defining elements of 4x4 matrix: each element is a function of x (more precisely: f(x^2)):

x       = sp.Symbol(x)
s1_14 = s1_21 = s1_34 = s1_41 = 0.0
s1_11 = A1*(x**2) + B1*x + C1
s1_12 = A2*(x**2) + B2*x + C2
s1_13 = A3*(x**2) + B3*x + C3
s1_22 = A1*(x**2) + B1*x + C1
s1_23 = A2*(x**2) + B2*x + C2
s1_24 = A3*(x**2) + B3*x + C3
s1_31 = A4*(x**2) + B4*x + C4
s1_32 = A5*(x**2) + B5*x + C5
s1_33 = A6*(x**2) + B6*x + C6
s1_42 = A4*(x**2) + B4*x + C4
s1_43 = A5*(x**2) + B5*x + C5
s1_44 = A6*(x**2) + B6*x + C6
s1    = sp.Matrix([[s1_11, s1_12, s1_13, s1_14], [s1_21, s1_22, s1_23, s1_24], [s1_31, s1_32, s1_33, s1_34], [s1_41, s1_42, s1_43, s1_44]])

# where A1,..., B1, ..., C1, ... are floats (computed in the same script before; all checked)
# for instance (they are small...):
# A1 = -8.832920000001262e-20
# A2 = 1.23393935e-16
# A3 = -1.841680000001036e-20
# A4 = 1.5494922020000018e-18
# A5 = 7.434205120000009e-19
# A6 = 5.154697260000002e-19
# B1 = -2.4092100000002057e-19
# B2 = -2.4352230000004432e-19
# B3 = -6.06127000000108e-20
# B4 = 3.555121236000002e-18
# B5 = 4.224810150000003e-18
# B6 = 1.2527727940000004e-18
# C1 = -1.5992580000002828e-19
# C2 = -1.656777000000101e-19
# C3 = -4.2327499999991735e-20
# C4 = 1.9217166420000004e-18
# C5 = 2.3215129980000005e-18
# C6 = 6.981400700000001e-19


# 2. calculate a determinant of that matrix and write it as a polynomial of x (the maximum order of that polynomial is then 8):

det = sp.simplify(s1.det())

# we can print the results; for coefficients A1...C6 above it is:
# det = (1.07479861449248e-87*x**10 + 8.00921385801962e-87*x**9 + 2.45757883595634e-86*x**8 + 3.97197994040761e-86*x**7 + 3.56434445101232e-86*x**6 + 1.68328171836074e-86*x**5 + 3.26853998044205e-87*x**4 + 4.42042971009037e-91*x**3 + 6.87690887243217e-92*x**2 - 1.67739461567147e-96*x - 2.66764671309489e-97)/(8.83292000000126e-20*x**2 + 2.40921000000021e-19*x + 1.59925800000028e-19)

# now for instance i do:

if det.is_rational_function():
    det_poly = sp.apart(det, x)
else:
    break

# what for the given coefficients gives:
# det_poly = 1.21681008601044e-68*x**8 + 5.74856653371868e-68*x**7 + 9.94041734367267e-68*x**6 + 7.44692074414804e-68*x**5 + 2.04344273909691e-68*x**4 + 2.1162503508103e-72*x**3 + 4.30019156172244e-73*x**2 - 7.97549194718549e-78*x + 1.0*(7.19921945884898e-101*x + 7.35103305249383e-101)/(8.83292000000126e-20*x**2 + 2.40921000000021e-19*x + 1.59925800000028e-19) - 1.66851240787895e-78


# combining sp.apart() with factor() or collect() does not change anything...

Aaron Meurer

unread,
Dec 22, 2017, 10:47:23 PM12/22/17
to sy...@googlegroups.com
I would consider this to be a bug in det(). It should try to compute
it in such a way that it returns a polynomial. This is somewhat
related to this issue https://github.com/sympy/sympy/issues/11082.

cancel() would normally be what you'd use to do this, but doesn't help
here either because it only cancels things if they cancel out exactly,
but because of the floating point coefficients, this does not happen.

apart() actually might work, if it didn't have its own issues
(https://github.com/sympy/sympy/issues/11795).

Can you open an issue for this in the issue tracker?

Unfortunately, I can't think of a good workaround for this right now.

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/ce698f69-2f11-40ae-a13c-b7194603ee3a%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

gosia olejniczak

unread,
Dec 28, 2017, 2:35:32 PM12/28/17
to sy...@googlegroups.com
Dear Aaron,

Thank you for your message and explanation!
I opened an issue:
https://github.com/sympy/sympy/issues/13805
and pasted your comments there;
best wishes,
gosia


> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/ce698f69-2f11-40ae-a13c-b7194603ee3a%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+unsubscribe@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.

Richard Fateman

unread,
Jan 8, 2018, 12:16:41 PM1/8/18
to sympy
Presumably  det() is using a method which requires exact division / cancellation
with polynomials, so it will work if all those floats are converted to rationals
and the arithmetic done exactly. The results will typically be numbers with
many many digits.

There are many papers on different methods for computing determinants
with symbolic entries. (e.g. fraction-free Gaussian elimination)
  There are at least 3 algorithms in Macsyma/ Maxima.

RJF

> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/ce698f69-2f11-40ae-a13c-b7194603ee3a%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.
Reply all
Reply to author
Forward
0 new messages