# possible bug in residue function

I compute the residue of a function in two ways, using the fact that it has simple poles. It seems residue method gives the wrong anser (psi) to be compared to correct answer (chi).

   br(x) = 1-x
q1,q2,q3,m = var('q1,q2,q3,m')
assume(q1>0,q2>0,q3>0,m>0)
q4 = (q1*q2*q3)^-1
k = 2
rho = [1,q4]
X = [var("x%d" % i) for i in range(k)]
chi1 = prod([ br(X[j]/m)/(br(X[j])*X[j]) for j in range(k)])
chi2 = prod([ prod([ br(q1*q2*X[i]/X[j])*br(q1*q3*X[i]/X[j])*br(q2*q3*X[i]/X[j]) for i in range(k) if i != j]) for j in range(k)])
chi3 = prod([ prod([ br(q1*X[i]/X[j])*br(q2*X[i]/X[j])*br(q3*X[i]/X[j])*br(q4*X[i]/X[j]) for i in range(k) if i != j]) for j in range (k)])
chi = (chi1*chi2/chi3).factor()
psi = chi
for xi,rhoi in zip(X,rho):
psi = psi.residue(xi==rhoi).combine().factor()
chi = (chi*(xi-rhoi)).factor().subs({xi: rhoi})


Now we can print chi and psi and see that they're different. Is this a bug of residue or series?

edit retag close merge delete

Sort by ยป oldest newest most voted

WorksForMe(TM) in Sage 9.1?beta0 ... when done manually:

sage: list(zip(X,rho))
[(x0, 1), (x1, 1/(q1*q2*q3))]
sage: psi.residue(x0==1).combine().factor()
(q1*q2*x1 - 1)*(q1*q3*x1 - 1)*(q2*q3*x1 - 1)*(q1*q2 - x1)*(q1*q3 - x1)*(q2*q3 - x1)*(m - x1)*(m - 1)*q1^2*q2^2*q3^2/((q1*q2*q3*x1 - 1)*(q1*q2*q3 - x1)*(q1*x1 - 1)*(q2*x1 - 1)*(q3*x1 - 1)*m^2*(q1 - x1)*(q2 - x1)*(q3 - x1)*(x1 - 1))
sage: (chi*(x0-1)).factor().subs(x0==1)
(q1*q2*x1 - 1)*(q1*q3*x1 - 1)*(q2*q3*x1 - 1)*(q1*q2 - x1)*(q1*q3 - x1)*(q2*q3 - x1)*(m - x1)*(m - 1)*q1^2*q2^2*q3^2/((q1*q2*q3*x1 - 1)*(q1*q2*q3 - x1)*(q1*x1 - 1)*(q2*x1 - 1)*(q3*x1 - 1)*m^2*(q1 - x1)*(q2 - x1)*(q3 - x1)*(x1 - 1))
sage: bool(psi.residue(x0==1).combine().factor()==(chi*(x0-1)).factor().subs(x0==
....: 1))
True
sage: psi.residue(x1==1/(q1*q2*q3)).combine().factor()
0
sage: (chi*(x1-1/(q1*q2*q3))).factor().subs(x1==1/(q1*q2*q3))
0
sage: bool(psi.residue(x1==1/(q1*q2*q3)).combine().factor()==(chi*(x1-1/(q1*q2*q3
....: ))).factor().subs(x1==1/(q1*q2*q3)))
True


I think that the problem resides in your two last lines of code:

    psi = psi.residue(xi==rhoi).combine().factor()
chi = (chi*(xi-rhoi)).factor().subs({xi: rhoi})


which scratch the previous values of psi and chi...

And, BTW, psi has a native definition in sage: from psi?:

Signature:      psi(x, *args, **kwds)
Docstring:
The digamma function, psi(x), is the logarithmic derivative of the
gamma function.

psi(x) = frac{d}{dx} log(Gamma(x)) =
frac{Gamma'(x)}{Gamma(x)}

We represent the n-th derivative of the digamma function with
psi(n, x) or psi(n, x).


Unless you have a pressing necessity of redefining it, it might be wiser to leave psi's definition untouched...

EDIT :

However, my last two lines are exactly as they shuold be: I want to iterate over the residue operation

Okay. Let's replace your last five lines by:

chi = [(chi1*chi2/chi3).factor()]
psi = chi.copy()
for xi,rhoi in zip(X,rho):
psi.insert(0,psi[0].residue(xi==rhoi).combine().factor())
chi.insert(0,(chi[0]*(xi-rhoi)).factor().subs({xi: rhoi}))


Then check each step:

sage: list(map(lambda u,v:bool(u==v), psi, chi))
[False, True, True]


The first iteration gives the expected result, the second cannot be proved True.

more

The equalities can be checked with:

list(map(lambda xi,rhoi:bool(psi.residue(xi==rhoi).combine().factor() == (chi*(xi-rhoi)).factor().subs({xi: rhoi})), X, rho))


HTH,

( 2020-01-22 06:35:27 +0200 )edit

Thanks for the tip about not redefining psi. However, my last two lines are exactly as they shuold be: I want to iterate over the residue operation. The problem I wanted to point out is that residue/series fails to give the correct result. Do you agree? Btw, do email notifications work for you?

( 2020-01-22 10:55:21 +0200 )edit

In particular, if you substitute e.g. q1=3, q2=5, q3=7, m=11 you should see that chi[0] and psi[0] are different.

( 2020-01-24 11:49:02 +0200 )edit

## Stats

Seen: 208 times

Last updated: Jan 22 '20