29
\$\begingroup\$

Challenge

Unlike the circumference of a circle (which is as simple as \$2\pi r\$), the circumference 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\$
  • 6
    \$\begingroup\$ I’m waiting for an answer that draws a bigger ellipse and counts the pixels \$\endgroup\$ – user Sep 29 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\$ – Mark Morgan Lloyd Sep 29 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\$ – Nick Sep 29 at 13:53
  • 3
    \$\begingroup\$ youtube.com/watch?v=5nW3nJhBHL0 \$\endgroup\$ – Alnitak Sep 29 at 16:01
  • 2
    \$\begingroup\$ @user competing for least serious answer codegolf.stackexchange.com/questions/211763/… \$\endgroup\$ – Stef Sep 30 at 14:30

24 Answers 24

13
\$\begingroup\$

Wolfram Language (Mathematica), 20 bytes

Perimeter[#~Disk~#]&

Try it online!

-2 bytes from @Roman (see comments)

| improve this answer | |
\$\endgroup\$
  • 5
    \$\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\$ – Michael Seifert Sep 30 at 13:47
  • 1
    \$\begingroup\$ -2 by removing N@ as allowed in the comments. \$\endgroup\$ – Roman Sep 30 at 19:04
12
\$\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.

| improve this answer | |
\$\endgroup\$
  • 1
    \$\begingroup\$ Nice solution! For the first one, it looks like k>>9or keeps enough precision. \$\endgroup\$ – xnor Sep 29 at 23:50
  • \$\begingroup\$ @xnor Indeed, thanks! \$\endgroup\$ – Anders Kaseorg Sep 30 at 0:22
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.

| improve this answer | |
\$\endgroup\$
  • 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 at 8:13
  • \$\begingroup\$ 23 bytes. Looks like 22 is possible... \$\endgroup\$ – Bubbler Sep 29 at 8:30
  • \$\begingroup\$ @Thanks again. This is going to take me some time to understand ;). \$\endgroup\$ – ovs Sep 29 at 8:45
8
\$\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)
| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ Can you please explain the "drawing" algorithm? \$\endgroup\$ – Stef Sep 30 at 16:10
  • \$\begingroup\$ @Stef done - I hope you approve! \$\endgroup\$ – Dominic van Essen Sep 30 at 16:19
  • \$\begingroup\$ Awesome! I thought this would require a Bresenham-like algorithm but this approach is even better. \$\endgroup\$ – Stef Sep 30 at 16:24
6
\$\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

| improve this answer | |
\$\endgroup\$
5
\$\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!

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ 57 bytes with a bit of rearrangement \$\endgroup\$ – Dominic van Essen Sep 29 at 8:29
  • \$\begingroup\$ @DominicvanEssen Thank you! -1 by rearranging a bit more. \$\endgroup\$ – Arnauld Sep 29 at 9:34
  • \$\begingroup\$ Nice. Very frustrating that R doesn't have +=... \$\endgroup\$ – Dominic van Essen Sep 29 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.

| improve this answer | |
\$\endgroup\$
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.

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ It could be 93 bytes as recursive function \$\endgroup\$ – Dominic van Essen Sep 29 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\$ – Dominic van Essen Sep 29 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 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\$ – Dominic van Essen Sep 29 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 at 16:29
4
\$\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\$.)

| improve this answer | |
\$\endgroup\$
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)
| improve this answer | |
\$\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.

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

| improve this answer | |
\$\endgroup\$
2
+150
\$\begingroup\$

APL (Dyalog Extended), 28 bytes

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

Try it online!

ovs's conversion to a train.

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

| improve this answer | |
\$\endgroup\$
2
\$\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.

| improve this answer | |
\$\endgroup\$
2
\$\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)
| improve this answer | |
\$\endgroup\$
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 ) } \$.

| improve this answer | |
\$\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!

| improve this answer | |
\$\endgroup\$
1
\$\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)$$

| improve this answer | |
\$\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.

| improve this answer | |
\$\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.

| improve this answer | |
\$\endgroup\$
  • \$\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\$ – Anders Kaseorg Oct 2 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 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 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\$ – Anders Kaseorg Oct 3 at 23:22
  • \$\begingroup\$ As for -MList::Util=sum, we count the command line flag as +17 bytes. \$\endgroup\$ – Anders Kaseorg Oct 3 at 23:32
1
\$\begingroup\$

CJam, 29 bytes

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

Try it online!

| improve this answer | |
\$\endgroup\$
1
\$\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
| improve this answer | |
\$\endgroup\$
0
\$\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.

| improve this answer | |
\$\endgroup\$
0
\$\begingroup\$

Japt, 35 bytes

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

Try it

| improve this answer | |
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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