Complete integral of the third kind using scipy/python and mathematica

I recently came to an issue which requires using the complete elliptical integral of the third kind. It was written as a single fuction in the form of PI(m,n). Honestly I am not familiar with this topic all I need is to calculate this value. I searched online and for python there are two packages available: scipy and mpmath. The statement of this calculation from these two packages are different: for mpmath (ellippi) it is simply "complete elliptical integral of the third kind" and for scipy (elliprj) "Symmetric elliptic integral of the third kind." These two statements already raises doubles for me. I tried mpmath first and somehow it produces complex numbers. And afterwards I tried scipy and it produces a tuple containing 4 numbers. And in the end I tried Mathematica (EllipticPi) and it only produces one value.

import scipy as scp
from mpmath import ellippi,ellipk
import numpy as np

r=1
rp=2
z=5
zp=7

gamma=zp-z
a=np.sqrt(gamma**2+(rp+r)**2)
c=np.sqrt(gamma**2+r**2)
n1=np.sqrt(2*r/np.abs(r-c))
k=np.sqrt(4*r*rp/a**2)
sn = ellippi(n1**2,k)

def ellippi(n, m):
  return scp.special.elliprf(0., 1. - m, 1.) + (n / 3.) * scp.special.elliprj(0., 1. - m, 1., 1. - n)

print(sn, ellippi(n1**2,k))

r = 1;
rp = 2;
z = 5;
zp = 7;
gamma = zp - z;
a = Sqrt[gamma^2 + (rp + r)^2];
b = rp + r;
c = Sqrt[gamma^2 + r^2];
n1 = Sqrt[2*r/Abs[(r - c)]];
k = Sqrt[4*r*rp/a^2];
N[EllipticPi[n1^2, k]]

My questions are:

1. what are the differences between mpmath.ellippi and scipy.special.elliprj
2. Should I use scipy pacakge and how should I transform the 4-number-tuple into one number?
3. Why scipy produce 4 numbers but Matematica only produce 1 number?

SciPy computes a more general integral (Carlson's). These integrals are ugly for general parameters, and, if all you need is the elliptic integral of third kind, you don't need to deal with these.

The relationship between what Mathematica gives you and Carlson's integrals (what scipy.special offers) is

def ellippi(n, m):
  return special.elliprf(                                                                                                                                                                                        
        0., 1. - m, 1.) + (n / 3.) * special.elliprj(0., 1. - m, 1., 1. - n)

which is also what mpmath gives called with two arguments, where n is the "characteristic", and m is the "parameter", i.e. the elliptic modulus squared.

Note that SciPy no longer supports Python 2.7, so running this code with older Python may result in 'scipy.special.ellipr*' modules not found.