33
\$\begingroup\$

Challenge

Unlike the circumference of a circle (which is as simple as \$2\pi r\$), 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.

enter image description here

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
\$\endgroup\$
16
  • 9
    \$\begingroup\$ I’m waiting for an answer that draws a bigger ellipse and counts the pixels \$\endgroup\$
    – user
    Sep 29, 2020 at 11:32
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Sep 29, 2020 at 13:51
  • 1
    \$\begingroup\$ @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 \$\endgroup\$ Sep 29, 2020 at 13:53
  • 3
    \$\begingroup\$ youtube.com/watch?v=5nW3nJhBHL0 \$\endgroup\$
    – Alnitak
    Sep 29, 2020 at 16:01
  • 2
    \$\begingroup\$ @user competing for least serious answer codegolf.stackexchange.com/questions/211763/… \$\endgroup\$
    – Stef
    Sep 30, 2020 at 14:30

27 Answers 27

18
\$\begingroup\$

Wolfram Language (Mathematica), 20 bytes

Perimeter[#~Disk~#]&

Try it online!

-2 bytes from @Roman (see comments)

\$\endgroup\$
2
  • 6
    \$\begingroup\$ 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. \$\endgroup\$ Sep 30, 2020 at 13:47
  • 2
    \$\begingroup\$ -2 by removing N@ as allowed in the comments. \$\endgroup\$
    – Roman
    Sep 30, 2020 at 19:04
15
\$\begingroup\$

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.

\$\endgroup\$
2
  • 1
    \$\begingroup\$ Nice solution! For the first one, it looks like k>>9or keeps enough precision. \$\endgroup\$
    – xnor
    Sep 29, 2020 at 23:50
  • \$\begingroup\$ @xnor Indeed, thanks! \$\endgroup\$ Sep 30, 2020 at 0:22
10
\$\begingroup\$

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)
\$\endgroup\$
3
  • \$\begingroup\$ Can you please explain the "drawing" algorithm? \$\endgroup\$
    – Stef
    Sep 30, 2020 at 16:10
  • \$\begingroup\$ @Stef done - I hope you approve! \$\endgroup\$ Sep 30, 2020 at 16:19
  • \$\begingroup\$ Awesome! I thought this would require a Bresenham-like algorithm but this approach is even better. \$\endgroup\$
    – Stef
    Sep 30, 2020 at 16:24
8
\$\begingroup\$

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.

\$\endgroup\$
3
  • 2
    \$\begingroup\$ Turns out using ⍳9 instead of ⍳99 is precise enough, and inlining x and using 1⊥ instead of +⌿∘ gives 25 bytes. \$\endgroup\$
    – Bubbler
    Sep 29, 2020 at 8:13
  • \$\begingroup\$ 23 bytes. Looks like 22 is possible... \$\endgroup\$
    – Bubbler
    Sep 29, 2020 at 8:30
  • \$\begingroup\$ @Thanks again. This is going to take me some time to understand ;). \$\endgroup\$
    – ovs
    Sep 29, 2020 at 8:45
7
\$\begingroup\$

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\$.)

\$\endgroup\$
0
7
\$\begingroup\$

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:

enter image description here

\$\endgroup\$
6
\$\begingroup\$

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!

\$\endgroup\$
3
  • \$\begingroup\$ 57 bytes with a bit of rearrangement \$\endgroup\$ Sep 29, 2020 at 8:29
  • \$\begingroup\$ @DominicvanEssen Thank you! -1 by rearranging a bit more. \$\endgroup\$
    – Arnauld
    Sep 29, 2020 at 9:34
  • \$\begingroup\$ Nice. Very frustrating that R doesn't have +=... \$\endgroup\$ Sep 29, 2020 at 9:55
4
\$\begingroup\$

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.

\$\endgroup\$
3
4
\$\begingroup\$

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.

\$\endgroup\$
5
  • \$\begingroup\$ It could be 93 bytes as recursive function \$\endgroup\$ Sep 29, 2020 at 12:08
  • \$\begingroup\$ ...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...) \$\endgroup\$ Sep 29, 2020 at 12:09
  • \$\begingroup\$ @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 \$\endgroup\$
    – Noodle9
    Sep 29, 2020 at 13:02
  • \$\begingroup\$ I think ceilingcat actually meant k?:2 and not k?k:2, so -2 bytes to 91 bytes... a new trick for me! \$\endgroup\$ Sep 29, 2020 at 15:57
  • \$\begingroup\$ @Dominic It is indeed k?:2 which is new to me too. Credited both of you for that byte - thanks! :-) \$\endgroup\$
    – Noodle9
    Sep 29, 2020 at 16:29
3
\$\begingroup\$

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)
\$\endgroup\$
3
\$\begingroup\$

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.

\$\endgroup\$
0
3
\$\begingroup\$

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)$$

\$\endgroup\$
3
\$\begingroup\$

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)
\$\endgroup\$
3
\$\begingroup\$

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
\$\endgroup\$
2
\$\begingroup\$

SageMath, 37 bytes

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

Try it online!

Uses the elliptic integral formulation.

\$\endgroup\$
2
\$\begingroup\$

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).

\$\endgroup\$
0
2
\$\begingroup\$

Charcoal, 52 bytes

≧×χφNθNηI×⁴ΣEEφE²∕⁺ιλφ₂⁺××θθ⁻Σι⊗₂Πι××ηη⁻⁻²Σι⊗₂⁻⊕ΠιΣι

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.

NθNη

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

I×⁴Σ

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.

EEφE²∕⁺ιλφ

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 ) } \$.

\$\endgroup\$
0
2
+150
\$\begingroup\$

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\$.

\$\endgroup\$
1
1
\$\begingroup\$

Symja, 35 bytes

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

Try It Online!

A port of the SageMath answer in Symja.

\$\endgroup\$
1
\$\begingroup\$

Japt, 35 bytes

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

Try it

\$\endgroup\$
1
\$\begingroup\$

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!

\$\endgroup\$
1
\$\begingroup\$

Pyth, 40 bytes

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

Try it online!

Just formula 5, like most other answers here.

\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
7
  • \$\begingroup\$ 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. \$\endgroup\$ Oct 2, 2020 at 9:06
  • \$\begingroup\$ 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. \$\endgroup\$
    – Kjetil S
    Oct 3, 2020 at 18:20
  • \$\begingroup\$ @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. \$\endgroup\$
    – Kjetil S
    Oct 3, 2020 at 18:24
  • \$\begingroup\$ 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. \$\endgroup\$ Oct 3, 2020 at 23:22
  • \$\begingroup\$ As for -MList::Util=sum, we count the command line flag as +17 bytes. \$\endgroup\$ Oct 3, 2020 at 23:32
1
\$\begingroup\$

CJam, 29 bytes

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

Try it online!

\$\endgroup\$
1
\$\begingroup\$

PARI/GP, 31 bytes

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

Attempt This Online!

\$\endgroup\$
1
\$\begingroup\$

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)\$.

\$\endgroup\$
1
  • \$\begingroup\$ elliprg also exists in scipy I notice. \$\endgroup\$
    – Simd
    Feb 19, 2023 at 4:49
0
\$\begingroup\$

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
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.