# Two polynomials multiply incorrectly only on my computer

I am using Docker Desktop on Windows 11 to run Sagemath with Macaulay package. The polynomials "F[0]" and "xd" below multiply into "dat" incorrectly (when compared to the correct result on another computer.)

I tried to obtain a minimal example illustrating my problem:

Source Code

import numpy as np
import sympy as sp
import itertools

alpha_2, alpha_01, x_0, x_1, x_2 = var('alpha_2 alpha_01 x_0 x_1 x_2')
y=np.array([x_0, x_1, x_2])

f=np.zeros(shape=(3),dtype=object)
f[0]=alpha_01*(x_0 + x_1 + x_2 )*(x_1 + x_2 ) + alpha_2*x_1*x_2 + (x_0*x_1 + x_0*x_2 + x_1*x_2 )*alpha_01 + 3*(x_0 + x_1 + x_2)^2
f[1]=alpha_01*(x_0 + x_1 + x_2 )*(x_0 + x_2 ) + alpha_2*x_0*x_2 + (x_0*x_1 + x_0*x_2 + x_1*x_2 )*alpha_01 + 3*(x_0 + x_1 + x_2)^2
f[2]=alpha_01*(x_0 + x_1 + x_2 )*(x_0 + x_1 ) + alpha_2*x_0*x_1 + (x_0*x_1 + x_0*x_2 + x_1*x_2 )*alpha_01 + 3*(x_0 + x_1 + x_2)^2

K=[4,1,1]
list2=[0,1,2]

F=np.zeros(shape=(3),dtype=object)
for p in range(3):
b=list2[p]
g=(f[b])^K[b]/y[b]^(2*K[b])*(1/K[b])
F[p]=(-g.subs({x_0:1})).simplify()    #F[i]=\hat{f}_i

A=sp.zeros(int(2),int(2))

for i, j in itertools.product(range(2),range(2)):
A[i,j]=diff( F[i+1], y[list2[j+1]] )

xd=det(A)

#https://doc.sagemath.org/html/en/reference/calculus/sage/calculus/test_sympy.html

print("F[0]=",F[0].expand())
print( " ---------------------------------- " )

print("xd=",xd.expand())

dat=F[0]*xd
dat=expand(dat)
print( " ---------------------------------- " )

print("DAT=F[0]*xd=",dat)


This program works fine on another computer. The result I get on my computer has a summand "-alpha_01^6x_1^11x_2^3". In a correct expansion of "dat", the summand "...+ -alpha_01^6x_1^11x_2^3+..." must not appear.

Note: I tried using polynomial rings instead of symbolic variables, but differentiation operation is not defined for polynomial rings.

edit retag close merge delete

1

It appears to be an issue with symbolic calculation. However, there seems to be no point to deal with symbolic expressions here, while you in fact are working with polynomial/rational functions. I therefore suggest to replace alpha_2, alpha_01, x_0, x_1, x_2 = var('alpha_2 alpha_01 x_0 x_1 x_2') with R.<alpha_2, alpha_01, x_0, x_1, x_2> = QQ[] and remove later simplify() and expand() calls as redundant. The issue should go away.

( 2022-04-12 23:49:04 +0200 )edit
1

Look at

sage: x, y = polygens(QQ, 'x,y')
sage: p = x * x + y
sage: p.derivative(x)
2*x

( 2022-04-13 09:03:12 +0200 )edit

As you recommended, I started with rational functions, and using derivatives and determinants that are already defined, I was able to obtain "dat" correctly! However, to obtain the coefficient of 1/x_1, I still need to import sympy and convert "dat" into sympy.core.mul.Mul, which makes me wonder if there is an easier way. To me, it seems like symbolic expressions are better in some respects and polynomial functions are better in other respects and it is possible to convert them into one another, e.g. via sympify().

( 2022-04-16 04:59:32 +0200 )edit

I've added an answer with more details,in particular explaining how to obtain the coefficient of 1/x_1. In fact, you don't need numpy and sympy here - Sage can handle everything more efficiently and in consistent manner.

( 2022-04-16 05:26:18 +0200 )edit

Sort by » oldest newest most voted

It appears to be an issue with symbolic calculation. However, there seems to be no point to deal with symbolic expressions here, while you in fact are working with polynomial/rational functions. It seems that it's more suitable to use LaurentPolynomialRing as a base class here. Also, it's better to leave dealing with matrices to Sage, not numpy. Here is an updated code:

R.<alpha_2, alpha_01, x_0, x_1, x_2> = LaurentPolynomialRing(QQ)
y=[x_0, x_1, x_2]

#f=np.zeros(shape=(3),dtype=object)
f=[0]*3
f[0]=alpha_01*(x_0 + x_1 + x_2 )*(x_1 + x_2 ) + alpha_2*x_1*x_2 + (x_0*x_1 + x_0*x_2 + x_1*x_2 )*alpha_01 + 3*(x_0 + x_1 + x_2)^2
f[1]=alpha_01*(x_0 + x_1 + x_2 )*(x_0 + x_2 ) + alpha_2*x_0*x_2 + (x_0*x_1 + x_0*x_2 + x_1*x_2 )*alpha_01 + 3*(x_0 + x_1 + x_2)^2
f[2]=alpha_01*(x_0 + x_1 + x_2 )*(x_0 + x_1 ) + alpha_2*x_0*x_1 + (x_0*x_1 + x_0*x_2 + x_1*x_2 )*alpha_01 + 3*(x_0 + x_1 + x_2)^2

K=[4,1,1]
list2=[0,1,2]

#F=np.zeros(shape=(3),dtype=object)
F=[0]*3
for p in range(3):
b=list2[p]
g=(f[b])^K[b]/y[b]^(2*K[b])*(1/K[b])
F[p]=(-g.subs({x_0:1}))    #F[i]=\hat{f}_i

A = Matrix(2,2,lambda i,j: diff(F[i+1],y[list2[j+1]]))

xd=det(A)

print("F[0]=",F[0])
print( " ---------------------------------- " )

print("xd=",xd)

dat=F[0]*xd
#dat=expand(dat)
print( " ---------------------------------- " )

print("DAT=F[0]*xd=",dat)


The coefficient of 1/x_1 (as a polynomial in other variables) in dat can be obtained as dat.coefficient(x_1^(-1)).

more

If an issue is confirmed, could you please open a ticket ?

( 2022-04-16 12:46:52 +0200 )edit

I do not see the described issue.

( 2022-04-16 18:26:59 +0200 )edit