c++cmath

C++ Elliptic Integral of first kind

I'm working on a program that makes calculations in C++; one of which involves using the complete elliptic integral of the first kind.

Python and Mathematica produce these results

Python:

from scipy import special
print(special.ellipk(0.1))
print(special.ellipk(0.2))
print(special.ellipk(0.3))

with output

1.6124413487202192
1.659623598610528
1.713889448178791

Mathematica

EllipticK[0.1]
EllipticK[0.2]
EllipticK[0.3]

same output

1.61244
1.65962
1.71389

But in C++

#include <stdio.h>
#include <cmath>

int main()
{   
    printf("0.1 %f \n", std::comp_ellint_1(0.1));
    printf("0.2 %f \n", std::comp_ellint_1(0.2));
    printf("0.3 %f \n", std::comp_ellint_1(0.3));
    return 0;
}

calculates different numbers

0.1 1.574746
0.2 1.586868
0.3 1.608049

It seems odd that there would be a mistake in C++'s standard cmath library. Is there a solution to this?

EDIT: When a small value (like 0.00001) is used, the numbers start to agree, so I think the issue is the calculation is an insufficient approximation. I am wondering if there is a version of this function for C++ that is as accurate as Mathematica or SciPy.

Solution

  • Elaborating on Steve's answer: there is no mistake or inaccuracy (select isn't broken). They are simply using different normalizations.

    Specifically, if you carefully read the documentation, C++'s std::comp_ellint_1(k) returns integral from 0 to pi/2 of dtheta/sqrt(1-k^2 sin^2 theta)

    whereas Python's scipy.special.ellipk(m) returns

    integral from 0 to pi/2 of dtheta/sqrt(1-m sin^2 theta)

    Note the k^2 versus m. So to get the result you want, you should call std::comp_ellint_1(std::sqrt(0.1)).

    Try it on godbolt

    This is consistent with your observation that the results agree more closely when the inputs are close to 0. If a number is very close to zero, then so is its square root, and both formulas yield a result very close to pi/2. If your input is really really small, then the difference will seem to vanish entirely due to underflow.

    This kind of discrepancy is unfortunately common in mathematics; there are no universally standardized definitions for mathematical functions, and it's not unusual to find variations. One that's always annoyed me: when defining the Fourier transform and inverse Fourier transform of an arbitrary function, there are several different options for where you could put the necessary factors of 2 pi. So if you look it up in four different textbooks, you may very well find it done in four different and incompatible ways. Thus you can't stop when you see the words "Fourier transform"; you have to look for the specific definition they are actually using, and adapt it if necessary to match what you want.