EllipticPi with complex arguments

[jbarth-ubhd]

Calling acb_elliptic_pi_inc(a_res, a_n, a_phi, a_m, quasiperiod_1, 52) like this (Mathematica syntax)

EllipticPi[+0.000000 -3.000000 Sqrt[-1], 
+0.769972 +0.680250 Sqrt[-1], 
-9.000000 -0.000000 Sqrt[-1]]  
= +0.575274 -0.173587 i

but Wolfram Mathematica states that

EllipticPi[+0.000000 - 3.000000 Sqrt[-1], 
+0.769972 +0.680250 Sqrt[-1],
 -9.000000 - 0.000000 Sqrt[-1]]
= 0.261115 + 0.768891 i

In my application (length of a cubic curve) the latter is much more reasonable, because Integral[b]-Integral[a] is non-complex.

PS: this issue is a private case and not associated with my profession.

[fredrik-johansson]

See also: https://github.com/fredrik-johansson/mpmath/issues/517

[fredrik-johansson]

Unfortunately, I have no idea how Mathematica defines this function with complex parameters.

[fredrik-johansson]

Can you provide some more information? What are the two values b, a such that Integral[b]-Integral[a] is non-complex?

[msalat]

I did want to know the length of a curve of type y = a x³ +b x² +c x +d

So I did with my raspberry pi

Integrate[Sqrt[1 + (3*a*x*x + 2*b*x + c)^2], x, 
 Assumptions -> {a \[Element] Reals, b \[Element] Reals, 
   c \[Element] Reals, x \[Element] Reals}]

This leads to a somewhat lengthy expression containing EllipticF(), EllipticE() and EllipticPi() (all with possibly complex arguments)

The result for a=4 and b=-4 and c=d=0 from x=0...1 should be around 1.63 (non-complex, estimated with numerical integration). With the function returned by Integrate[] it does so in mathematica.

My C program to calculate the mathematica generated expression is:

#include <cstdio>
#include <cmath>
#include <complex>
#include <assert.h>
#include "arb.h"
#include "acb_elliptic.h"

// N[Integrate[Sqrt[1 + (12*x*x - 8*x)^2], {x, 0, 1}]]
// mit a=4 und b=-4 von 0..1: 1.63218

/* Integrate[Sqrt[1 + (3*a*x*x + 2*b*x + c)^2], x, 
 Assumptions -> {a \[Element] Reals, b \[Element] Reals, 
   c \[Element] Reals, x \[Element] Reals}]
*/

using namespace std;

#define Sqrt sqrt
#define Power pow
#define Complex complex<double>
#define ArcSin asin

const int quasiperiod_1=0;

complex<double> EllipticF(complex<double>phi, complex<double>m) {
    acb_t a_m, a_phi, a_res;
    acb_init(a_m); acb_init(a_phi), acb_init(a_res);
    acb_set_d_d(a_phi, phi.real(), phi.imag());
    acb_set_d_d(a_m  ,   m.real(),   m.imag());
    acb_elliptic_f(a_res, a_phi, a_m, quasiperiod_1, 52);
    arf_struct *real=arb_midref(acb_realref(a_res));
    arf_struct *imag=arb_midref(acb_imagref(a_res));
    complex<double> res(arf_get_d(real, ARF_RND_NEAR), arf_get_d(imag, ARF_RND_NEAR));
    acb_clear(a_m);
    acb_clear(a_phi);
    acb_clear(a_res);
    printf("EllipticF[%+lf %+lf Sqrt[-1], %+lf %+lf Sqrt[-1]] = %+lf %+lf i\n", phi.real(), phi.imag(), m.real(), m.imag(), res.real(), res.imag());
    return res;
}

complex<double> EllipticE(complex<double>phi, complex<double>m) {
    acb_t a_m, a_phi, a_res;
    acb_init(a_m); acb_init(a_phi), acb_init(a_res);
    acb_set_d_d(a_phi, phi.real(), phi.imag());
    acb_set_d_d(a_m  ,   m.real(),   m.imag());
    acb_elliptic_e_inc(a_res, a_phi, a_m, quasiperiod_1, 52);
    arf_struct *real=arb_midref(acb_realref(a_res));
    arf_struct *imag=arb_midref(acb_imagref(a_res));
    complex<double> res(arf_get_d(real, ARF_RND_NEAR), arf_get_d(imag, ARF_RND_NEAR));
    acb_clear(a_m);
    acb_clear(a_phi);
    acb_clear(a_res);
    printf("EllipticE[%+lf %+lf Sqrt[-1], %+lf %+lf Sqrt[-1]] = %+lf %+lf i\n", phi.real(), phi.imag(), m.real(), m.imag(), res.real(), res.imag());
    return res;
}

complex<double> EllipticPi(complex<double>n, complex<double>phi, complex<double>m) {
    acb_t a_n, a_m, a_phi, a_res;
    acb_init(a_n);
    acb_init(a_phi);
    acb_init(a_m);
    acb_init(a_res);
    acb_set_d_d(a_n  ,   n.real(),   n.imag());
    acb_set_d_d(a_phi, phi.real(), phi.imag());
    acb_set_d_d(a_m  ,   m.real(),   m.imag());
    acb_elliptic_pi_inc(a_res, a_n, a_phi, a_m, quasiperiod_1, 52);
    // das hilft auch nicht: acb_printn(a_res, 50, 0); flint_printf("\n");
    arf_struct *real=arb_midref(acb_realref(a_res));
    arf_struct *imag=arb_midref(acb_imagref(a_res));
    complex<double> res(arf_get_d(real, ARF_RND_NEAR), arf_get_d(imag, ARF_RND_NEAR));
    acb_clear(a_n);
    acb_clear(a_phi);
    acb_clear(a_m);
    acb_clear(a_res);
    printf("EllipticPi[%+lf %+lf Sqrt[-1], %+lf %+lf Sqrt[-1], %+lf %+lf Sqrt[-1]] = %+lf %+lf i\n", n.real(), n.imag(), phi.real(), phi.imag(), m.real(), m.imag(), res.real(), res.imag());
    /*
OOPS diese Funktion hier sagt:
EllipticPi[+0.000000 -3.000000 Sqrt[-1], +0.769972 +0.680250 Sqrt[-1], -9.000000 -0.000000 Sqrt[-1]] 
= +0.575274 -0.173587 i -- Fehler auch in arb-2.17.0 UND auch in python-mpmath

Wolfram Mathematica sagt:
EllipticPi[+0.000000 - 3.000000 Sqrt[-1], +0.769972 + 
  0.680250 Sqrt[-1], -9.000000 - 0.000000 Sqrt[-1]]
= 0.261115 + 0.768891 I
*/

    return res;
}

double integral_len_cubic(double x, double a, double b, double c) { // d unnötig

complex<double>q1=Sqrt(Complex(0,3)*a + Power(b,2) - 3*a*c);
complex<double>q2=Sqrt(Power(b,2) - 3*a*(Complex(0,1) + c));
complex<double>q3=Sqrt(((q1 - q2)*(b + q2 + 3*a*x))/((q1 + q2)*(b - q2 + 3*a*x)));
complex<double>q4=Sqrt((q2*(-b + q1 - 3*a*x))/((q1 + q2)*(-b + q2 - 3*a*x)));
complex<double>q5=Sqrt(-((q2*(b + q1 + 3*a*x))/((-q1 + q2)*(-b + q2 - 3*a*x))));
complex<double>q6=EllipticF(ArcSin(q3),Power(q1 + q2,2)/Power(q1 - q2,2));
complex<double>q7=EllipticPi((q1 + q2)/(q1 - q2),ArcSin(q3), Power(q1 + q2,2)/Power(q1 - q2,2));
complex<double>q8=EllipticE(ArcSin(q3),Power(q1 + q2,2)/Power(q1 - q2,2));
double q9=1 + Power(c,2) + 4*Power(b,2)*Power(x,2) + 12*a*b*Power(x,3) + 9*Power(a,2)*Power(x,4) + 2*c*x*(2*b + 3*a*x);
complex<double>q10=Power(b - q2 + 3*a*x,2);
complex<double>q11=q10*(q1 + q2)*q3*q4*q5*q6;
complex<double>q12=-2.*(-Power(b,2) + 2.*b*q2 - q1*q2)*q6 + 8.*b*q2*q7 + Power(q1 - q2,2)*q8;
complex<double>q13=q10*q12*(q1 + q2)*q3*q4*q5;
complex<double>q14=q13 + q2*(-q1 + q2)*(b - q1 + 3*a*x)*(b + q1 + 3*a*x)*(b + q2 + 3*a*x);
complex<double>q15=q10*(q1 + q2)*q3*q4*q5;
complex<double>q16=q12*q15 + q2*(-q1 + q2)*(b - q1 + 3*a*x)*(b + q1 + 3*a*x)*(b + q2 + 3*a*x);

complex<double> r= (2.*(-(Power(b,2)*q16) + 3.*a*c*q16 + 18.*Power(a,2)*q15*q6 - 
        6.*a*Power(b,2)*c*q15*q6 + 18.*Power(a,2)*Power(c,2)*q15*q6 + 
        4.*Power(b,3)*q15*((b - q2)*q6 + 2.*q2*q7) - 
        12.*a*b*c*q15*((b - q2)*q6 + 2.*q2*q7)))/
    (81.*Power(a,3)*(-q1 + q2)*Sqrt((Power(b,2) - 3.*a*(Complex(0,1) + c))*q9)) + 
   (Sqrt(q9)*(b + 3*a*x))/(9.*a);
    printf("%+lf %+lfi\n", r.real(), r.imag());
   return r.real();
}

int main(int argc, char **argv) {
    flint_printf("Computed with arb-%s\n", arb_version);
    complex<double> r=integral_len_cubic(1, 4, -4, 0)-integral_len_cubic(0, 4, -4, 0);
   printf("%lf +%lf i\n", r.real(), r.imag());
}

[msalat]

PS: as far as I remember EllipticE+EllipticF do not differ from mathematica output (at least not for the input arguments in this context)

[fredrik-johansson]

Mathematica may be producing a symbolic integral tailored to its own pecular definition of EllipticPi.

The difficulties might go away if you express your integral in terms of Carlson forms instead of the Legendre elliptic integrals.

[fredrik-johansson]

Consider this algorithm: https://www.ams.org/journals/mcom/2002-71-237/S0025-5718-01-01333-3/

[laolux]

Just for reference, Mathematica output seems to be version dependent. For 11.2 I get the same output as Wolfram Alpha, i.e.

EllipticPi[+0.000000 - 3.000000 Sqrt[-1], 
+0.769972 +0.680250 Sqrt[-1], 
-9.000000 - 0.000000 Sqrt[-1]]
= 0.261115 + 0.768891 i

So it might not be a good idea to try to reproduce Mathematica's result.