# Circumference of an ellipse

## Challenge

Unlike the circumference of a circle (which is as simple as $$\2\pi r/extract_tex]), the circumference (arc length) of an ellipse is hard. Given the semi-major axis $$\a\$$ and semi-minor axis $$\b\$$ of an ellipse (see the image below, from Wikipedia), calculate its circumference. By definition, you can assume $$\0 < b \le a\$$ for input values. The output value must be within $$\10^{-6}\$$ relative error from the expected answer for the given test cases. Standard rules apply. The shortest code in bytes wins. ## Formulas Relevant information can be found on Wikipedia and MathWorld. $$\C\$$ is the value of the circumference; $$\e\$$ and $$\h\$$ are helper values. The last two are Ramanujan's approximations, the first of which (the crossed-out one) does not meet the error requirements. The second approximation formula (Equation 5) barely does (verification) for up to $$\a=5b\$$ (which is also the upper limit of the test cases, so you can use it for your answer). \require{enclose} \\ \begin{align} e &= \sqrt{1-\frac{b^2}{a^2}} \\ C &= 4aE(e) = 4a\int^{\pi/2}_{0}{\sqrt{1-e^2 \sin^2 \theta} \;d\theta} \tag{1} \\ C &= 2 \pi a \left(1-\sum^{\infty}_{n=1}{\left(\frac{(2n-1)!!}{(2n)!!}\right)^2 \frac{e^{2n}}{2n-1}}\right) \tag{2} \\ h &= \frac{(a-b)^2}{(a+b)^2} \\ C &= \pi (a + b) \left( 1 + \sum^{\infty}_{n=1} { \left( \frac{(2n-1)!!}{2^n n!} \right)^2 \frac{h^n}{(2n-1)^2} } \right) \tag{3} \\ C &= \pi (a + b) \sum^{\infty}_{n=0} { \binom{1/2}{n}^2 h^n } \tag{4} \\ \enclose{horizontalstrike}{C} &\enclose{horizontalstrike}{\approx \pi \left( 3(a+b) - \sqrt{(3a+b)(a+3b)} \right)} \\ C &\approx \pi (a+b) \left( 1+ \frac{3h}{10 + \sqrt{4-3h}} \right) \tag{5} \end{align} ## Test cases All the values for C (circumference) are calculated using Equation 4 with 1000 terms, and presented with 10 significant figures. a b C 1 1 6.283185307 1.2 1 6.925791195 1.5 1 7.932719795 2 1 9.688448220 3 1 13.36489322 5 1 21.01004454 20 10 96.88448220 123 45 556.6359936  • I’m waiting for an answer that draws a bigger ellipse and counts the pixels – user Sep 29, 2020 at 11:32 • Please can I highlight that a number of the proposed solutions use characters that aren't in the 7-bit ASCII character set, so it's not accurate to count their expression in "bytes": it should be characters or codepoints, some of which require several bytes for their composition. Sep 29, 2020 at 13:51 • @MarkMorganLloyd If you're talking about languages like APL (Dyalog Unicode), you may be interested in this meta post, you'll find that they frequently use a special character set and that number of bytes == number of characters Sep 29, 2020 at 13:53 • youtube.com/watch?v=5nW3nJhBHL0 Sep 29, 2020 at 16:01 • @user competing for least serious answer codegolf.stackexchange.com/questions/211763/… – Stef Sep 30, 2020 at 14:30 ## 27 Answers # Wolfram Language (Mathematica), 20 bytes Perimeter[#~Disk~#]&  Try it online! -2 bytes from @Roman (see comments) • It's amusing that even though EllipticE and Ellipsoid are built-ins, this is still the most economical way to do it. Also, that's a clever way to get around the first argument being the central coordinates. Sep 30, 2020 at 13:47 • -2 by removing N@ as allowed in the comments. Sep 30, 2020 at 19:04 # 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. • Nice solution! For the first one, it looks like k>>9or keeps enough precision. – xnor Sep 29, 2020 at 23:50 • @xnor Indeed, thanks! Sep 30, 2020 at 0:22 # 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)  • Can you please explain the "drawing" algorithm? – Stef Sep 30, 2020 at 16:10 • @Stef done - I hope you approve! Sep 30, 2020 at 16:19 • Awesome! I thought this would require a Bresenham-like algorithm but this approach is even better. – Stef Sep 30, 2020 at 16:24 # APL (Dyalog Unicode), 2825 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. • Turns out using ⍳9 instead of ⍳99 is precise enough, and inlining x and using 1⊥ instead of +⌿∘ gives 25 bytes. Sep 29, 2020 at 8:13 • 23 bytes. Looks like 22 is possible... Sep 29, 2020 at 8:30 • @Thanks again. This is going to take me some time to understand ;). – ovs Sep 29, 2020 at 8:45 # 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\$$.) # x87 machine code, 6559 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: # 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! • 57 bytes with a bit of rearrangement Sep 29, 2020 at 8:29 • @DominicvanEssen Thank you! -1 by rearranging a bit more. Sep 29, 2020 at 9:34 • Nice. Very frustrating that R doesn't have +=... Sep 29, 2020 at 9:55 # 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. • quick 1 byte savings: Try it online! Sep 29, 2020 at 12:53 • @Jonah Ah, so obvious now! Thanks! Sep 29, 2020 at 12:56 • @Jonah or alternatively [:o.1#.+*i.@9(2^~^~*0.5!~[)-%+ Sep 29, 2020 at 13:05 # 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. • It could be 93 bytes as recursive function Sep 29, 2020 at 12:08 • ...but I don't know if it's safe to assume that a non-specified function argument always initializes to zero (it seems to do so on TIO...) Sep 29, 2020 at 12:09 • @DominicvanEssen GCC is initializing \k\ to zero for each invocation of \f\ in the same program run so I believe we're ok - thanks! :D Sep 29, 2020 at 13:02 • I think ceilingcat actually meant k?:2 and not k?k:2, so -2 bytes to 91 bytes... a new trick for me! Sep 29, 2020 at 15:57 • @Dominic It is indeed k?:2 which is new to me too. Credited both of you for that byte - thanks! :-) Sep 29, 2020 at 16:29 # 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)  # 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. # 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)$$ # MATL, 19 bytes y/U_Q.5t_hlbZh*YPE*  ### 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$/extract_tex] is the eccentricity, \[ e = \sqrt{1 - b^2/a^2},$ and $$\E$/extract_tex] 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)  # APL(NARS), 21 chars a←⎕⋄{√+/×⍨a×1 2○⍵}∫○2  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π). It is possible calculate more digit in that integral but solution is more long, as in  ⎕FPC←3000⋄a←⎕⋄300⍕{√+/×⍨a×1 2○⍵}∫○2v ⎕: 1 1 6.283185307179586476925286766559005768394338798750211641949889184615632812572417997256069650684234135964296173026564613 294187689219101164463450718816256962234900568205403877042211119289245897909860763928857621951331866892256951296467 573566330542403818291297133846920697220908653296426787214520498282547 300⍕○2v 6.283185307179586476925286766559005768394338798750211641949889184615632812572417997256069650684234135964296173026564613 294187689219101164463450718816256962234900568205403877042211119289245897909860763928857621951331866892256951296467 573566330542403818291297133846920697220908653296426787214520498282547  # SageMath, 37 bytes lambda a,b:4*a*elliptic_ec(1-b*b/a/a)  Try it online! Uses the elliptic integral formulation. # 05AB1E, 2221 20 bytes ÆnIOn/3*D4s-tT+/>IOžqP  Implements the fifth formula. Input as a pair $$\[a,b]\$$. -1 byte thanks to @ovs. 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, $$\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). # 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 ) } \$$. # APL (Dyalog Extended), 28 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\$$. # Symja, 35 bytes f=N(4*#1*EllipticE(1-#2*#2/#1/#1))&  Try It Online! A port of the SageMath answer in Symja. # Japt, 35 bytes MP*ºH=3*(U-V ²/(U±V)/(A+(4-H/U ¬ +U  Try it # 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! # Pyth, 40 bytes A,hQeQJc^-GH2^+GH2**.n0+GHhc*3J+T@-4*3J2  Try it online! Just formula 5, like most other answers here. # Perl 5, 70 bytes sub{mys;maps+=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, 7465 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. • Our usual rules require you to include use List::Util 'sum'; in the byte count, no? That said, map sqrt((_[0]*cos_)**2+(_[1]*sin_)**2) saves 4 bytes. Oct 2, 2020 at 9:06 • Maybe? I swapped the first answer to my four byte longer answer. I thought that since List::Util has been a core module for many years (since 2002), all normal Perl installations have it. Seems a matter of "administration" for a language which functions are true built-ins and which are semi-built-in in core modules that doesnt require extra installations. I've seen C answers without #include <math.h> and Javascript answers regularly omits f= in function definitions. Oct 3, 2020 at 18:20 • @AndersKaseorg Thx for the cos-sin-tip. And it can even loose five extra bytes with map sqrt+(_[0]*cos)**2+(_[1]*sin)**2. I get a warning about sqrt not having parens. Sin and cos uses _ without args. Oct 3, 2020 at 18:24 • There’s no exemption for matters of “administration”. A C program can still run (with a warning) without #include <math.h>, and we allow function expressions without f= (unless the function is recursive). Your program does not run without use List::Util 'sum';. The issue is not whether it exists to be imported; the issue is that must be explicitly imported to be used even when it exists. Oct 3, 2020 at 23:22 • As for -MList::Util=sum, we count the command line flag as +17 bytes. Oct 3, 2020 at 23:32 # CJam, 29 bytes {_:+_P*@:-@d/_*3*_4\-mqA+/)*}  Try it online! # PARI/GP, 31 bytes f(a,b)=4*a*ellE((1-b^2/a^2)^.5)  Attempt This Online! # 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)\$$. • elliprg also exists in scipy I notice. – Simd Feb 19, 2023 at 4:49 # 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