APL(NARS), 18 chars

a←⎕⋄{∣<a×1 2○⍵}∫○2

{∣<a×1 2○⍵}∫○2 would mean the integral of the function {∣<a×1 2○⍵} in the interval 0..2π where ○2 return 2π.

I find one unducumentated instruction < that given input array of 2 numbers a b, return the number complex aJb and so can apply the function norm ∣ to <a×1 2○⍵ (it means 3 chars less). test:

      a←⎕⋄{∣<a×1 2○⍵}∫○2
⎕:
      1 1
6.283185307 
      a←⎕⋄{∣<a×1 2○⍵}∫○2
⎕:
      123 45
556.6359936 

it is based on the parametric representation of ellipse something as

phi(t)=(a cos(t),b sin(t)) t∊[0,2π) and the formula for find the lenght of the curve $$L(\varphi)=\int_0^{2\pi} \left| \varphi'(t) \right|dt $$

where phi'(t)=(-a sin(t), b cos(t) ) t∊[0,2π).

So we would have |phi'(t)|=sqrt(a^2 sin(t)^2+ b^2 cos(t)^2)

Pheraps it is possible calculate more digits in that integral but solution is more long, as in

  ⎕FPC←300×4⋄a←⎕⋄300⍕{∣<a×1 2○⍵}∫○2v
⎕:
  1 1
6.283185307179586476925286766559005768394338798750211641949889184615632812
  572417997256069650684234135964296173026564613294187689219101164463450
  718816256962234900568205403877042211119289245897909860763928857621951
  331866892256951296467573566330542403818291297133846920697220908653296
  426787214520498282547
  ⎕FPC←300×4⋄a←⎕⋄300⍕{∣<a×1 2○⍵}∫○2v
⎕:
  123 45
556.6359936189455567443781051948463024126034035502138146906597334028357954
  725986409879935974903682461041935521442813245838254811371393443438747
  420629698312290554610067139047841633148575060990440709560230306713908
  847002927362101581753119688993061698742679070849585365528914366426158
  76601356181848049458166

note the 2v at end of the line.

Python/mpmath, 51 bytes

lambda a,b:8*elliprg(0,a*a,b*b)
from mpmath import*

Attempt This Online!

Using Carlson's symmetric integrals – as seen here – not only removes the need to sort input, it saves 1 byte off the bog-standard \$4aE(e)\$ formula. As given in DLMF 19.30.5, \$L=8R_G(0,a^2,b^2)\$.

PARI/GP, 31 bytes

f(a,b)=4*a*ellE((1-b^2/a^2)^.5)

Attempt This Online!

APL (Dyalog Extended), ²⁸ 26 bytes

○+×1+∘(⊢÷10+∘√4-⊢)3×2*⍨-÷+

Try it online!

ovs's conversion to a train.

-2 using √.

APL (Dyalog Extended), 35 bytes

{h←3×2*⍨⍺(-÷+)⍵⋄(○⍺+⍵)×1+h÷10+√4-h}

Try it online!

Uses Equation 4.

Longer than the other APL answer because there's more than one usage of \$h\$.

Arn, 22 bytes

┴þ5‡Ô縄”R¤ËíÜç›WðÙÝÁ*

Try it! A pretty good approximation, but not exact for the larger values. Uses the crossed out formula (which I assume was removed due to the innacuracy). For any wondering, I managed to get the non-crossed out formula 5 to 33 bytes, but I couldn't figure out how to shorten it (and it was even less accurate than this one).

Explained

Unpacked: pi*(3*(+\)-:/(*3+:})*+3*:}

pi                     Variable; first 20 digits of π
  *
    (
          3
        *
          (+\)         Folded sum ([a, b] -> a + b)
      -
        :/             Square root
            (
                  _    Variable; initialized to STDIN; implied
                *
                  3
              +
                  _    Implied
                :}     Tail
            )
          *
              _        Implied
            +
                3
              *
                  _    Implied
                :}
                       Ending parentheses implied

CJam, 29 bytes

{_:+_P*@:-@d/_*3*_4\-mqA+/)*}

Try it online!

x87 machine code, 65 59 53 bytes

00000000: d9c1 d9c1 dec1 d9ca dee9 d8c8 d9c1 d8c8  ................
00000010: def9 6a03 8bf4 de0c ff04 df04 d9c1 dee9  ..j.............
00000020: d9fa 8304 06de 04de f9d9 e8de c1d9 ebde  ................
00000030: c9de c95e c3                             ...^.

Listing:

D9 C1       FLD   ST(1)             ; load a to ST
D9 C1       FLD   ST(1)             ; load b to ST
DE C1       FADD                    ; a + b
D9 CA       FXCH  ST(2)             ; save result for end 
DE E9       FSUB                    ; a - b 
D8 C8       FMUL  ST(0), ST(0)      ; ST ^ 2 
D9 C1       FLD   ST(1)             ; copy a + b result to ST 
D8 C8       FMUL  ST(0), ST(0)      ; ST ^ 2 
DE F9       FDIV                    ; calculate h 
6A 03       PUSH  3                 ; load const 3 
8B F4       MOV   SI, SP            ; SI to top of CPU stack
DE 0C       FIMUL WORD PTR[SI]      ; ST = h * 3 
FF 04       INC   WORD PTR[SI]      ; 4 = 3 + 1 
DF 04       FILD  WORD PTR[SI]      ; load const 4 
D9 C1       FLD   ST(1)             ; load 3h to ST
DE E9       FSUB                    ; 4 - 3h 
D9 FA       FSQRT                   ; sqrt(ST) 
83 04 06    ADD   WORD PTR[SI], 6   ; 10 = 4 + 6 
DE 04       FIADD WORD PTR[SI]      ; ST + 10 
DE F9       FDIV                    ; 3h / ST 
D9 E8       FLD1                    ; load const 1 
DE C1       FADD                    ; ST + 1 
D9 EB       FLDPI                   ; load PI 
DE C9       FMUL                    ; * PI 
DE C9       FMUL                    ; * ( a + b ) from earlier 
5E          POP   SI                ; restore CPU stack 
C3          RET                     ; return to caller 

Callable function, input a and b in ST(0) and ST(1). Output in ST(0). Implements Ramanujan's 2nd approximation (eq 5) in full hardware 80-bit extended precision.

Test program:

Charcoal, 52 bytes

≧×χφＮθＮηＩ×⁴ΣＥＥφＥ²∕⁺ιλφ₂⁺××θθ⁻Σι⊗₂Πι××ηη⁻⁻²Σι⊗₂⁻⊕ΠιΣι

Try it online! Link is to verbose version of code. Works by approximating the line integral for a quadrant. The default precision is unfortunately only ~5 significant figures so the first four bytes are needed to increase the precision to ~7 significant figures. Further increases are possible for the same byte count but then it becomes too slow to demonstrate on TIO. Explanation:

≧×χφ

Increase the number of pieces \$ n \$ in which to divide the quadrant from \$ 1,000 \$ to \$ 10,000 \$. ≧×φφ would increase it to \$ 1,000,000 \$ but that's too slow for TIO.

ＮθＮη

Input the ellipse's axes \$ a \$ and \$ b \$.

Ｉ×⁴Σ

After calculating the approximate arc length of each piece into which the quadrant was subdivided, take the sum, multiply by \$ 4 \$ for the whole ellipse and output the result.

ＥＥφＥ²∕⁺ιλφ

Create a list of pieces of the quadrant. In the ellipse equation \$ \left ( \frac x a \right ) ^ 2 + \left ( \frac y b \right ) ^ 2 = 1 \$ we can set \$ \left ( \frac {x_i} a \right ) ^ 2 = \frac i n \$ and \$ \left ( \frac {y_i} b \right ) ^ 2 = 1 - \frac i n \$. Given a piece index \$ i \$ we want to calculate the distance between \$ ( x_i, y_i ) \$ and \$ ( x _{i+1}, y_{i+1} ) \$. For each \$ i \$ we calculate \$ j = \frac i n \$ and \$ k = \frac {i+1} n \$ and loop over the list.

₂⁺××θθ⁻Σι⊗₂Πι××ηη⁻⁻²Σι⊗₂⁻⊕ΠιΣι

The distance \$ \sqrt { ( a \sqrt k - a \sqrt j ) ^ 2 + ( b \sqrt { 1 - j } - b \sqrt { 1 - k } ) ^ 2 } \$ expands to \$ \sqrt { a^2 \left ( j + k - 2 \sqrt { j k } \right ) + b^2 \left ( (1 - j) + (1 - k) - 2 \sqrt { (1 - j) (1 - k) } \right ) } \$ which expands to \$ \sqrt { a^2 \left ( j + k - 2 \sqrt { j k } \right ) + b^2 \left ( 2 - (j + k) - 2 \sqrt { 1 + j k - (j + k) } \right ) } \$.

MATL, 19 bytes

y/U_Q.5t_hlbZh*YPE*

Try it online! Or verify all test cases.

Formula used

This is based on formula (1) from the challenge description, \[ C = 4a\int^{\pi/2}_{0}{\sqrt{1-e^2 \sin^2 \theta} ;d\theta} = 4 a\,E(e), \] where \$e\$ is the eccentricity, \[ e = \sqrt{1 - b^2/a^2}, \] and \$E\$ is the complete elliptic integral of the second kind. This integral can be expressed in terms of Gauss' hypergeometric function, \${}_2F_1\$, as follows: \[ E(e) = \tfrac{\pi}{2} \;{}_2F_1 \left(\tfrac12, -\tfrac12; 1; e^2 \right). \] Combining the above gives the formula used in the code: \[ C = 2\pi a \;{}_2F_1 \left(\tfrac12, -\tfrac12; 1; 1 - b^2/a^2 \right). \]

Code explanation

y       % Implicit inputs: a, b. Duplicate from below
        % STACK: a, b, a
/       % Divide
        % STACK: a, b/a
U_Q     % Square, negate, add 1
        % STACK: a, 1-(b/a)^2
.5t_h   % Push 0.5, duplicate, negate, concatenate
        % STACK: a, 1-(b/a)^2, [0.5, -0.5]
1       % Push 1
        % STACK: a, 1-(b/a)^2, [0.5, -0.5], 1
b       % Bubble up in the stack
        % STACK: a, [0.5, -0.5], 1, 1-(b/a)^2
Zh      % Hypergeometric function, 2F1
        % STACK: a, 2F1([0.5, -0.5], 1, 1-(b/a)^2)
*       % Multiply
        % STACK: a * 2F1([0.5, -0.5], 1, 1-(b/a)^2)
YPE     % Push pi, multiply by 2
        % STACK: a * 2F1([0.5, -0.5], 1, 1-(b/a)^2), 2*pi
*       % Multiply. Implicit display
        % STACK: 2*pi*a * 2F1([0.5, -0.5], 1, 1-(b/a)^2)

Perl 5, 70 bytes

sub{my$s;map$s+=sqrt+($_[0]*cos)**2+($_[1]*sin)**2,0..1570795;4e-6*$s}

Try it online!

Perl 5, 78 bytes

sub f{($a,$b)=@_;$H=3*(($a-$b)/($a+=$b))**2;3.141593*$a*(1+$H/(10+sqrt 4-$H))}

With the a+=b trick stolen from the Javascript answer.

Try it online!

Or this one which is 13 bytes less (but uses core module List::Util)

Perl 5 -MList::Util=sum, 74 65 65+16 bytes

sub f{4e-6*sum map sqrt+($_[0]*cos)**2+($_[1]*sin)**2,0..1570795}

Try it online!

Which numerically calculates a variant of formula (1).

I was surprised this worked with sin and cos of integers up to 1570795 ≈ 500000π. But the tests in the question in "Try it online" has relative error < 0.000001. Guess sin²(the integers) is "averaged out" good enough.

Pyth, 40 bytes

A,hQeQJc^-GH2^+GH2**.n0+GHhc*3J+T@-4*3J2

Try it online!

Just formula 5, like most other answers here.

R, 60 57 bytes

function(a,b,c=a+b,h=3*(a-b)^2/c)pi*(c+h/(10+(4-h/c)^.5))

Try it online!

Straightforward implementation of Ramanujan's 2nd approximation (eq 5).

Rather sadly, this approximation comes out as much more concise than a more-interesting different approach prompted by the comments: 'draw' a big ellipse, and measure around the edge of it (unfortunately counting the actual pixels wasn't going to work...):

R, 90 65 62 bytes

Edit: -3 bytes by calculating hypotenuse length using abs value of complex number

function(a,b,n=1e5)sum(4*abs(diff(b*(1-(0:n/n)^2)^.5)+1i*a/n))

Try it online!

How? (ungolfed code):

circumference_of_ellipse=
function(a,b                # a,b = axes of ellipse
n=1e6){                     # n = number of pixels to 'draw' across 'a' axis
x=a*0:n/n                   # x coordinates = n pixels from 0 to a
y=b*(1-(x/a)^2)^.5)         # y coordinates = to satisfy (x/a)^2 + (y/b)^2 =1
                            # we could actually draw the (quarter) ellipse here
                            # with 'plot(x,y)'
step_y=diff(y)              # step_y = change in y for each step of x
step_x=a/n                  # step_x = size of each step of x
h=(step_y^2+step_x^2)^.5    # h=hypotenuse of triangle formed by step_y & step_x
sum(4*h)                    # sum all the hypotenuses and multiply by 4
                            # (since we only 'drew' a quarter of the ellipse)

Haskell, 73 bytes

e a b=(a+b)*pi*(1+3*l/(10+sqrt(4-3*l))+3*l^5/2^17)where l=((a-b)/(a+b))^2

Experimenting with an improved version of (5):

$$E(a,b) = \pi (a+b) \left( 1 + \frac{3h^2}{10 + \sqrt{4-3h^2}} + \frac{3h^{10}}{2^{17}}\right)$$

Wolfram Language (Mathematica), 20 bytes

Perimeter[#~Disk~#]&

Try it online!

-2 bytes from @Roman (see comments)

Wolfram Language (Mathematica), 25 24 bytes

4EllipticE[1-(#2/#)^2]#&

Try it online!

-1 thanks to @AndersKaseorg

Note that Mathematica uses a different convention for elliptic integrals, hence the square root disappears.

Python 3, 68 67 bytes

f=lambda a,b,k=2:k>>9or(1-b*b/a/a)*(k-4+3/k)/k*f(a,b,k+2)+6.28319*a

Try it online!

An exact infinite series, given sufficiently accurate values of \$2\pi \approx 6.28319\$ and \$\infty \approx 9\$.

69 68 bytes

f=lambda a,b,k=0:k//7*.785398*a*(8-k)or f(a+b,2*(a*b)**.5,k*b/a/2+4)

Try it online!

Another exact series, given sufficiently accurate values of \$\frac\pi4 \approx .785398\$ and \$8 \approx 7\$. This one converges extremely quickly, using just five recursive calls for each test case! The recursion exactly preserves the invariant value

$$\left(1 + \frac{kb}{8a}\right)C(a, b) - \frac{kb}{8a}C(a + b, 2\sqrt{a b}),$$

which can then be approximated as \$(1 - \frac k8)2\pi a\$ when \$a, b\$ become sufficiently close.

Ruby, 63 bytes

->a,b,h=1r*(a-b)/a+=b{3.141593*a*((154+53*h*=h)*h*h/1e4+h/4+1)}

Try it online!

A direct port of @Arnauld's JavaScript answer is shorter (58 bytes). However, I like the 63-byter above because it differs from other approaches in that it's a cubic polynomial: no square roots, no infinite series.

This excellent review lists nearly 40 different methods for approximating the circumference of an ellipse, with graphs of the relative error in each approximation as a function of \$b/a\$. Inspection of the graphs shows that only a few of the listed methods are capable of satisfying the required tolerance of \$10^{-6}\$ for all test cases. Since several answers here had already explored 'Ramanujan II' (eq. (5)), I decided to look at the Padé approximations 'Padé 3/2' and 'Padé 3/3'.

A Padé approximant is a rational function with coefficients chosen so as to match the largest possible number of terms in a known power series. In this case, the relevant power series is the infinite sum that appears in eq. (4). The Padé 3/2 and Padé 3/3 approximants for this series are mathematically straightforward (see the review linked above) but not suited to code golf. Instead, an approximation to the approximants is obtained by least-squares fitting. The resulting cubic polynomial (with truncated coefficients), as implemented in the code, is $$ 0.0053h^3 + 0.0154h^2+0.25h+1. $$ Note that this function is overfitted to the test cases, partly because of the truncation and partly because the fit was optimised using only those values of \$h=(a-b)^2/(a+b)^2\$ that occur in the test cases. (Consequently, Math::PI cannot be substituted in place of 3.141593, despite having the same byte count, without yielding relative errors above the \$10^{-6}\$ threshold for the two test cases for which \$b/a=1/2\$.)

Jelly, 20 bytes

I÷S²3×÷ạ4½+⁵Ʋ$‘×S×ØP

A monadic Link accepting a pair of [a, b] which yields the result of formula 5.

Try it online!

I thought formula 4 would be the way to go, but only got 21:

9Ḷ.c×⁹I÷S*⁸¤²ʋ€×ØP×SS

Try it online!

C (gcc), ^{97 92} 91 bytes

Saved 4 5 bytes thanks to Dominic van Essen!!!
Saved 2 bytes thanks to ceilingcat!!!

float f(a,b,k)float a,b,k;{k=k?:2;k=k>999?1:(1-b*b/a/a)*(k-4+3/k)/k*f(a,b,k+2)+6.283185*a;}

Try it online!

Port of Anders Kaseorg's Python answer.

Japt, 35 bytes

MP*ºH=3*(U-V ²/(U±V)/(A+(4-H/U ¬ +U

Try it

J, 31 30 bytes

-1 byte thanks to Jonah!

[:o.1#.+*i.@9*:@(^~*0.5!~[)-%+

Try it online!

Essentially a J port of @ovs's APL solution.

APL (Dyalog Unicode), 28 25 23 bytes

Thanks to Bubbler for -5 bytes!

Assumes ⎕IO←0.

f←○1⊥+×9(×⍨*×.5!⍨⊢)∘⍳⍨-÷+

Try it online!

This calculates

$$ \pi \cdot \sum_{n=0}^{8} (a+b) \cdot \left( h^{\prime n} \binom{1/2}{n} \right) ^2 \qquad h^\prime = {{a-b}\over{a+b}} $$

which is a good enough approximation using the 4th formula. For the explanation the function will be split into two. f is the main function and g calculates \$ \left( \alpha^{\prime n} \binom{1/2}{n} \right) ^2 \$ for \$n\$ from \$0\$ to \$\omega-1\$:

g ← (×⍨*×.5!⍨⊢)∘⍳
f ← ○1⊥+×9g⍨-÷+

Starting with a f b from the right:

-÷+ calculates \$h^\prime = (a-b)÷(a+b)\$.
g⍨ is g commuted => 9 g⍨ h' ≡ h' g 9. g returns a vector of the 9 values of \$\left( h^{\prime n} \binom{1/2}{n} \right) ^2\$.
+× multiplies \$a + b\$ to this vector.
1⊥ converts the resulting vector from base 1, which is the same as summing the vector.
○ multiplies the resulting number by \$\pi\$.

Now to h' g 9:

⍳ is an index generator, with ⎕IO←0, ⍳9 results in the vector 0 1 ... 8.
The remaining train ×⍨*×.5!⍨⊢ is now called with \$h^\prime\$ as a left argument and the vector \$v = (0,1, \cdots, 8)\$ as a right argument:

.5!⍨⊢ is the commuted binomial coefficient called with the vector v on its right and \$0.5\$ on its left. This calculates \$\binom{1/2}{n}\$ for all \$n \in v\$.
*× multiplies this vector element-wise with \$h^\prime * n\$ (\$*\$ denotes exponentiation).
×⍨ is commuted multiplication, which given only a right argument, seems to use this as left and right argument? and squares the vector element-wise.

JavaScript (ES7),  59  56 bytes

Saved 2 bytes thanks to @DominicvanEssen

a=>b=>Math.PI*((h=3*(a-b)**2/(a+=b))/(10+(4-h/a)**.5)+a)

Try it online!

MathGolf, 20 bytes

-ëΣ_¬/²3*_4,√♂+/)π**

Port of my 05AB1E answer, and thus also implements a modification of the fifth formula.

Try it online.

Explanation:

-                     # b-a
 ëΣ                   # a+b
   _                  # Duplicate
    ¬                 # Rotate stack: b-a,a+b,a+b → a+b,b-a,a+b
     /                # Divide
      ²               # Square
       3*             # Multiply by 3
         _            # Duplicate
          4,          # Subtract from 4
            √         # Square-root
             ♂+       # Add 10
               /      # Divide
                )     # Increment by 1
                 π*   # Multiply by PI
                   *  # Multiply by the a+b we've duplicated
                      # (after which the entire stack is output implicitly as result)

05AB1E, 22 21 20 bytes

ÆnIOn/3*D4s-tT+/>IOžqP

Implements the fifth formula. Input as a pair \$[a,b]\$.

-1 byte thanks to @ovs.

Try it online or verify all test cases.

Explanation:

Æ                     # Reduce the (implicit) input-pair by subtraction: a-b
 IO                   # Push the input-pair again and sum it: a+b
   /                  # Divide them by one another: (a-b)/(a+b)
    n                 # Square it: ((a-b)/(a+b))²
     3*               # Multiply it by 3: ((a-b)/(a+b))²*3
       D              # Duplicate that
        4α            # Take the absolute difference with 4: |((a-b)/(a+b))²*3-4|
          t           # Take the square-root of that: sqrt(|((a-b)/(a+b))²*3-4|)
           T+         # Add 10: sqrt(|((a-b)/(a+b))²*3-4|)+10
             /        # Divide the duplicate by this:
                      #  (a-b)²/(a+b)²*3/(sqrt(|((a-b)/(a+b))²*3-4|)+10)
              >       # Increase it by 1:
                      #  (a-b)²/(a+b)²*3/(sqrt(|((a-b)/(a+b))²*3-4|)+10)+1
               IO     # Push the input-sum again: a+b
                 žq   # Push PI: 3.141592653589793
                   P  # Take the product of the three values on the stack:
                      #  ((a-b)²/(a+b)²*3/(sqrt(|((a-b)/(a+b))²*3-4|)+10)+1)*(a+b)*π
                      # (after which the result is output implicitly)

Note that I use \$\left|3h-4\right|\$ instead of \$4-3h\$ in my formula to save a byte, but given the constraints \$0<b\leq a\$, \$h\$ will be: \$0\leq h<1\$, and thus \$3h\$ will be at most \$2.999\dots\$.
I also use \$h=\left(\frac{a-b}{a+b}\right)^2\$ instead of \$h=\frac{(a-b)^2}{(a+b)^2}\$ to save another byte (thanks to @ovs).

Symja, 35 bytes

f=N(4*#1*EllipticE(1-#2*#2/#1/#1))&

Try It Online!

A port of the SageMath answer in Symja.

SageMath, 37 bytes

lambda a,b:4*a*elliptic_ec(1-b*b/a/a)

Try it online!

Uses the elliptic integral formulation.