Complete Elliptic Integral of the First Kind

in C99
An implementation of
Arithmetic-Geometric Mean
$a_{n+1} &=& \frac{1}{2}(a_n + b_n)\\ b_{n+1} &=& \sqrt{a_n b_n}$

Jan. 8 2023
Takayuki HOSODA

The arithmetic-geometric mean can be used to compute - among others - logarithms, complete and incomplete elliptic integrals of the first and second kind, and Jacobi elliptic functions.

agm.c

NAME
    _agm -- the Arithmetic-Geometric Mean

LIBRARY
    Math library (libm, -lm)\n\

SYNOPSIS
    #include <math.h>
    #include <complex.h>
    #include <float.h>

    double complex
    _agm(double complex a, double complex b);

RETURN VALUE
    _agm returns the arithmetic-geometric mean agm(a, b) of two numbers a and b.

#include <math.h>
#include <complex.h>
#include <float.h>

double complex _agm(double complex a, double complex b);

/** The arithmetic-geometric mean
 _agm(a, b) returns the arithmetic-geometric mean agm(a, b) of two numbers a and b.
 a_{n+1} &=& \frac{1}{2}(a_n + b_n)\\
 b_{n+1} &=&  \sqrt{a_n b_n}
 */
double complex
_agm(a, b)
    double complex a;
    double complex b;
{
    double         e;
    double complex m;
    if ((cabs(a) == 0) || (cabs(b) == 0) || (a == -b)) return 0;
    e = 2 * DBL_EPSILON * fmax(cabs(a), cabs(b));
    do {
        m = csqrt(a * b);
        a = (a + b) * 0.5;
        b = ((cabs(a + m) == cabs(a - m)) && (cimag(m / a) > 0)) || (cabs(a + m) < cabs(a - m)) ? -m : m;
    } while (cabs(a - b) > e);
    return (a + b) * 0.5;
}

#ifdef DEBUG
#include <stdio.h>
#include <stdlib.h>
char *help = "Name\n\
    agm -- returns the arithmetic-geometric mean agm(a, b) of two numbers a and b.\n\
    Rev. 1.21 (Jan. 8, 2023) (c) 2022, Takayuki HOSODA, Finetune co., ltd.\n\
Usage\n\
    agm(a, b) with real arguments : agm a b\n\
    agm(a + bi, c + di) with complex arguments : agm a b c d\n\
    agm(1, z) for pm3D plot: agm range_re range_im samples\n";
int
main(int argc, char *argv[]) {
    double complex a, b, m;
    double         x, y;
    int            n, r, j;
    if(argc == 5) {
        a = atof(argv[1]) + I * atof(argv[2]);
        b = atof(argv[3]) + I * atof(argv[4]);
        m = _agm(a,  b); 
        printf("agm(%.16g %+.16gi, %.16g %+.16gi) ~= %.16g %+.16gi\n",
            creal(a), cimag(a), creal(b), cimag(b), creal(m), cimag(m));
    } else if (argc == 3) {
        a = atof(argv[1]);
        b = atof(argv[2]);
        m = _agm(a, b); 
        if (cimag(m) == 0) {
            printf("agm(%.16g, %.16g) ~= %.16g\n", creal(a), creal(b), creal(m));
        } else {
            printf("agm(%.16g, %.16g) ~= %.16g%+.16gi\n", creal(a), creal(b), creal(m), cimag(m));
        }
    } else if (argc == 4) {
        x = atof(argv[1]);
        y = atof(argv[2]);
        n = atof(argv[3]);
        if (n < 10) n = 10;
        printf("#agm(1, z) computed by \"agm Rev.1.21 (Jan. 8, 2023) (c) 2022, Takayuki HOSODA\".\n");
        for (r = 0 ; r <= n; r++) {
            for (j = 0 ; j <= n; j++) {
                a = 1.0;
                b = (-x * r + (n - r) * x) / (double)n + I * ((-y * j + (n - j) * y) / (double)n);
                m = _agm(a, b);
                printf("%.16g\t%.16g\t%.16g\t%.16g\n",creal(b), cimag(b), cabs(m), carg(m));
            }
            printf("\n");
        }
    } else {
        printf("%s\n", help);
    }
    return(0);
}
#endif

Calculation example of _agm()

agm( 1       ,  2     ) ~=  1.456791031046907
agm(-2   +2 i,  1     ) ~=  0.0126420344317752 +1.323914959272823i
agm(-3   +2 i, -2 +2 i) ~= -2.48482321791346 +2.012328388475829i
agm( 1   +2 i, -3 -4 i) ~= -2.051752595659329 +0.07245869868588142i
agm(-1.5 -2 i,  3 +4 i) ~=  1.956981679536812 -0.14596685837510370i

Calculation results by WolframAlpha

agm(1, 2)               ~= 1.4567910310469068691864323832650819749738639432213055907941723832...
agm(-2 + 2 i, 1)        ~= 0.012642034431775224871144882985098306376405625767985282751139254... 
                         + 1.3239149592728229110180317045523087521571829295502353536455198... i
agm(-3 + 2 i, -2 + 2 i) ~= -2.48482321791345966540888302608455218648014162219542779589616405...
                         + 2.01232838847582894914613621798033481637378492177117070798341113... i
agm(1 + 2 i, -3 - 4 i)  ~= -2.051752595659329572693583864239458955398878779423072275607602...
                          + 0.07245869868588149813029197687616500840652284108968969453100464... i
agm(-1.5 - 2i, 3 + 4 i) ~= 1.95698167953681205998137584021754054093691534541639368664817292...
                         - 0.145966858375103771971664986222696913809981322179346508693130731... i



SEE ALSO

Arithmetic-Geometric Mean — for HP 42S
Complete Elliptic Integral of the First Kind — in C99
Conductor-backed Coplanar Waveguide
Calculation of π by using arithmetic geometric mean


REFERENCE

Complete Elliptic Integral of the First Kind — Wolfram MathWorld
Arithmetic-Geometric Mean — Wolfram MathWorld
Louis. V King, "ON THE DIRECT NUMERICAL CALCULATION OF ELLIPTIC FUNCTIONS AND INTEGRALS", Cambridge at the university press, 1924









© 2000 Takayuki HOSODA.