AGM Series Hole 1: Calculate the Arithmetic–Geometric Mean

Question

This question was inspired by this HNQ.

About the series

This question is now part of a series about the AGM method. This first post in the series will be about actually calculating the AGM. You may treat this like any other code golf challenge, and answer it without worrying about the series at all. However, there is a leaderboard across all challenges.

What is the Arithmetic–Geometric Mean

The Arithmetic–Geometric Mean of two numbers is defined as the number that repeatedly taking the arithmetic and geometric means converges to. Your task is to find this number after some n iterations.

Clarifications

• You take three numbers, a, b, n in any reasonable format.
• For n iterations, take the arithmetic and geometric mean of a and b and set those to a and b.
• For two numbers a and b, the arithmetic mean is defined as (a + b) / 2.
• The geometric mean is defined as √(a * b).
• a and b should be approaching each other.
• Then, output both a and b.
• You don't have to worry about float imprecision and such.
• This is code-golf so shortest code in bytes wins!

Test Cases

[0, [24, 6]] -> [24, 6]    
[1, [24, 6]] -> [15.0, 12.0]
[2, [24, 6]] -> [13.5, 13.416407864998739]
[5, [24, 6]] -> [13.458171481725616, 13.458171481725616]
[10, [100, 50]] -> [72.83955155234534, 72.83955155234534]

The next one is 1/Gauss's Constant:
[10, [1, 1.41421356237]] -> [1.198140234734168, 1.1981402347341683]

Leaderboard

Stolen from Martin's series.

The following snippet will generate a leaderboard across all challenges of the series.

To make sure that your answers show up, please start every answer with a headline, using the following Markdown template:

# Language Name, N bytes

where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:

# Ruby, <s>104</s> <s>101</s> 96 bytes

/* Configuration */

var QUESTION_IDs = [66068]; // Obtain this from the url
// It will be like http://XYZ.stackexchange.com/questions/QUESTION_ID/... on any question page
var ANSWER_FILTER = "!.FjwQBrX2KXuFkv6p2lChi_RjzM19";

/* App */

var answers = [], page = 1, currentQ = -1;

function answersUrl(index) {
  return "http://api.stackexchange.com/2.2/questions/" +  QUESTION_IDs.join(";") + "/answers?page=" + index + "&pagesize=100&order=desc&sort=creation&site=codegolf&filter=" + ANSWER_FILTER;
}


function getAnswers() {
  $.ajax({
    url: answersUrl(page++),
    method: "get",
    dataType: "jsonp",
    crossDomain: true,
    success: function (data) {
      answers.push.apply(answers, data.items);
      if (data.has_more) getAnswers();
      else process();
    }
  });
}

getAnswers();

var SIZE_REG = /\d+(?=[^\d&]*(?:<(?:s>((?!>).)*<\/s>|((?!>).)+>)[^\d&]*)*$)/;
var NUMBER_REG = /\d+/;
var LANGUAGE_REG = /^#*\s*([^\n,]+)(?=,)/;//

function shouldHaveHeading(a) {
  var pass = false;
  var lines = a.body_markdown.split("\n");
  try {
    pass |= /^#/.test(a.body_markdown);
    pass |= ["-", "="]
              .indexOf(lines[1][0]) > -1;
    pass &= LANGUAGE_REG.test(a.body_markdown);
  } catch (ex) {}
  return pass;
}

function shouldHaveScore(a) {
  var pass = false;
  try {
    pass |= SIZE_REG.test(a.body_markdown.split("\n")[0]);
  } catch (ex) {}
  if (!pass) console.log(a);
  return pass;
}

function getAuthorName(a) {
  return a.owner.display_name;
}

function getAuthorId(a) {
  return a.owner.user_id;
}

function process() {
  answers = answers.filter(shouldHaveScore)
                   .filter(shouldHaveHeading);
  answers.sort(function (a, b) {
    var aB = +(a.body_markdown.split("\n")[0].match(SIZE_REG) || [Infinity])[0],
        bB = +(b.body_markdown.split("\n")[0].match(SIZE_REG) || [Infinity])[0];
    return aB - bB
  });

  var users = {};
  answers.forEach(function (a) {
    var headline = a.body_markdown.split("\n")[0];
    var question = QUESTION_IDs.indexOf(a.question_id);
    var size = parseInt((headline.match(SIZE_REG)||[0])[0]);
    var language = headline.match(LANGUAGE_REG)[1];
    var user = getAuthorName(a);
    var userId = getAuthorId(a);
    if (!users[userId]) users[userId] = {name: user, nAnswer: 0, answers: []};
    if (!users[userId].answers[question]) {
      users[userId].answers[question] = {size: Infinity};
      users[userId].nAnswer++;
    }
    if (users[userId].answers[question].size > size) {
      users[userId].answers[question] = {size: size, link: a.share_link}
    }
  });
  
  
  var sortedUsers = [];
  for (var userId in users)
    if (users.hasOwnProperty(userId)) {
      var user = users[userId];
      user.score = 0;
      user.completedAll = true;
      for (var i = 0; i < QUESTION_IDs.length; ++i) {
        if (user.answers[i])
          user.score += user.answers[i].size;
        else
          user.completedAll = false;
      }
      sortedUsers.push(user);
    }  
  
  sortedUsers.sort(function (a, b) {
    if (a.nAnswer > b.nAnswer) return -1;
    if (b.nAnswer > a.nAnswer) return 1;
    return a.score - b.score;
  });
  
  var place = 1;
  for (var i = 0; i < sortedUsers.length; ++i) {
    var user = sortedUsers[i];
    var row = '<tr><td>'+ place++ +'.</td><td>'+user.name+'</td>';
    for (var j = 0; j < QUESTION_IDs.length; ++j) {
      var answer = user.answers[j];
      if (answer)
        row += '<td><a href="'+answer.link+'">'+answer.size+'</a></td>';
      else
        row += '<td class="missing"></td>';
    }
    row += '<td></td>';
    if (user.completedAll)
      row += '<td class="total">'+user.score+'</td>';
    else
      row += '<td class="total missing">'+user.score+'</td>';
    row += '</tr>';
    $("#users").append(row);
  }
}

body { text-align: left !important}

#leaderboard {
  width: 500px; 
}

#answer-list {
  padding: 10px;
  width: 290px;
  float: left;
}

#language-list {
  padding: 10px;
  width: 290px;
  float: left;
}

table thead {
  font-weight: bold;
}

table td {
  padding: 5px;
}

td.total {
  font-weight: bold;
  text-align: right;
}

td.missing {
  background: #bbbbbb;
}

<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script>
<link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b">
<div id="leaderboard">
  <h2>Leaderboard</h2>
  <p>
    Missing scores are shown as grey cells. A grey total indicates that the user has not participated in all challenges and is not eligible for the overall victory yet.
  </p>
  <table class="_user-list">
    <thead>
      <tr><td></td><td>User</td>
        <td><a href="http://codegolf.stackexchange.com/q/66068">#1</a></td>
        <td></td><td>Total</td>
      </tr>
    </thead>
    <tbody id="users">

    </tbody>
  </table>
</div>
<table style="display: none">
  <tbody id="answer-template">
    <tr><td>{{NAME}}</td><td>{{LANGUAGE}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr>
  </tbody>
</table>
<table style="display: none">
  <tbody id="language-template">
    <tr><td>{{LANGUAGE}}</td><td>{{NAME}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr>
  </tbody>
</table>

Answer 1

CJam, 16 bytes

{{_2$+2/@@*mq}*}

Takes input on the stack as a b n where a and b are doubles. Online demo

Answer 2

Dyalog APL, 22 21 15 bytes

.5∘(+.×,×.*⍨)⍣⎕

Takes (a,b) as right argument, and prompts for n:

(

  +.× dot product of 0.5 and the right argument

 , followed by

  ×.*⍨ "dot power" of the right argument and 0.5^*

)⍣⎕ applied numeric-prompt times.

^* "dot power" is like dot product, but using multiplication and power instead of plus and multiplication, as follows:

      _n
A ×.*⍨ B is ∏ B_i^A = ∏ B₁^AB₂^{A
      i=1}

-3 bytes thanks to ngn.

---

Old version:

{((+/÷≢),.5*⍨×/)⍣⍺⊢⍵}

Takes n as left argument and a b as right argument.

⊢⍵ On the RightArg
(...)⍣⍺ recalculate LeftArg times
(+/÷≢) sum divided by tally
, followed by
.5*⍨×/ the square root of the product.

All the test cases:

      f←{((.5×+/),.5*⍨×/)⍣⍺⊢⍵}
      0 1 2 5 10 10 f¨ (24 6)(24 6)(24 6)(24 6)(100 50)(1,2*.5)
┌────┬─────┬────────────────┬───────────────────────┬───────────────────────┬───────────────────────┐
│24 6│15 12│13.5 13.41640786│13.45817148 13.45817148│72.83955155 72.83955155│1.198140235 1.198140235│
└────┴─────┴────────────────┴───────────────────────┴───────────────────────┴───────────────────────┘

Answer 3

TI-BASIC, 22 bytes

Input N
For(I,1,N
{mean(Ans),√(prod(Ans
End
Ans

Does exactly what the algorithm says. Takes N from the prompt, and A and B through Ans as a two-element list.

If N is 0, the For( loop is skipped entirely.

Answer 4

JavaScript ES7, 48 43 bytes

-5 thanks to Downgoat!

f=(n,a,b)=>n?f(n-1,(a+b)/2,(a*b)**.5):[a,b]

Very simple recursive function.

Answer 5

MATLAB/Octave, 69 65 bytes

function [a,b]=r(a,b,n)
for i=1:n;j=(a+b)/2;b=(a*b)^.5;a=j;end

Answer 6

Jelly, non-competing

9 bytes This answer is non-competing, since it uses features that postdate the challenge.

SH;P½¥ðṛ¡

Try it online!

How it works

SH;P½¥ðṛ¡    Input: x (vector) -- y (repetitions)

SH           Take the sum (S) of x and halve (H) the result.
   P½        Take the product (P) of x and the square root (½) of the result.
     ¥       Combine the last two instructions in a dyadic chain.
  ;          Concatenate the results to the left and to the right.
      ð      Push the preceding, variadic chain; begin a new, dyadic chain.
       ṛ     Return the right argument (y).
        ¡    Repeat the pushed chain y times.

Answer 7

C++, 108 102 100 bytes

Thank you to @RetoKoradi and @AlexA for saving me 6 bytes.

This is non-competitive, because C++ is not a good golfing language. Did this for fun :)

#include<cmath>
std::string f(float a,float b,int n){return n==0?a+" "+b:f((a+b)/2,sqrt(a*b),n-1);}

This is a simple recursion function, very similar to the JS answer.

Answer 8

Seriously, 11 bytes

,p`;π√@æk`n

Hex Dump:

2c70603be3fb40916b606e

Try it online

Explanation:

,                    Read in the list as [n,a,b]
 p                   pop list to yield: n [a,b]
  `      `n          Push a quoted function and run it n times.
   ;                 Duplicate [a,b] pair
    π√               Compute its product and square root it (GM)
      @              Swap the other copy of the pair to the top
       æ             Compute its mean.
        k            Compile the stack back into a list.

Answer 9

K5, 15 bytes

Very literal:

{(+/x%2;%*/x)}/

In action:

 {(+/x%2;%*/x)}/[0; 24 6]
24 6
 {(+/x%2;%*/x)}/[5; 24 6]
1.345817e1 1.345817e1

Unfortunately, this does not work in oK because that interpreter does not currently support projection (currying) of adverbs. Works in the real k5.

In oK, it would currently be necessary to wrap the definition in a lambda:

  {x{(+/x%2;%*/x)}/y}[5; 24 6]
13.4582 13.4582

Answer 10

J, 18 13 bytes

-:@+/,%:@*/^:

Usage:

   agm =: -:@+/,%:@*/^:
   5 agm 24 6
13.4582 13.4582

Answer 11

Japt, 24 bytes 25 33

Saved 9 7 bytes thank to @ETHproductions

Uo r@[VW]=[V+W /2(V*W q]

Takes advantage of ES6 destructuring.

Try it online

Ungolfed && Explanation

Uo r@[VW]=[V+W /2(V*W q]

       // Implicit: U: 1st input, V: 2nd input, W: 3rd input
Uo     // Range from 0 to 1st input
r@     // Loop over range
  [V,W]=    // Set 2nd and 3rd input to...
   [V+W /2,   // Add 2nd and 3rd inputs, divide by 2
   (V*W q]    // Multiple 2nd and 3rd inputs, find square root
            // Set's to the above respectively 
       // Implicit: return [V,W]

Answer 12

Pyth, 12

u,.OG@*FG2EQ

Test Suite

Explanation

u,.OG@*FG2EQ    ##  implicit: Q = eval(input())
u         EQ    ##  reduce eval(input()) times, starting with Q
                ##  the reduce lambda has G as the previous value and H as the next
  .OG           ##  arithmetic mean of last pair
     @*FG2      ##  geometric mean of last pair, uses *F to get the product of the list
                ##  and @...2 to get the square root of that
 ,              ##  join the two means into a two element list

Answer 13

Matlab, 54 bytes

function x=f(x,n)
for k=1:n
x=[mean(x) prod(x)^.5];end

Example:

>> f([24 6], 2)
ans =
  13.500000000000000  13.416407864998739

Answer 14

Minkolang v0.14, 23 bytes

Try it here!

$n[$d+2$:r*1Mi2%?!r]$N.
$n                      C get all input C
  [                ]    C pop N; repeat inner N times C
   $d                   C duplicate stack [1,2] => [1,2,1,2] C
     +                  C add top two elements C
      2$:               C divide by two C
         r              C reverse stack (get the other two) C
          *             C multiply them together C
           1M           C take square root C
             i2%?!r     C reverse the stack if an odd step number C
                    $N  C output stack
           1M           C take square root C
             i          C get step in for loop C

Answer 15

Python 3, 65 55 bytes

Thanks to mathmandan for pointing out a shorter version using the lambda operator.

f=lambda a,b,n:f((a+b)/2,(a*b)**.5,n-1)if n else(a,b)

My original version:

def f(a,b,n):
 if n:f((a+b)/2,(a*b)**.5,n-1)
 else:print(a,b)

To my chagrin, a recursive function (a la the JavaScript and C++ answers) was shorter than a simple for loop.

Answer 16

Pyth, 15 bytes

u,^*FG.5csG2vzQ

Answer 17

Mathematica, 31 30 bytes

Saved one byte thanks to Martin Büttner.

{+##/2,(1##)^.5}&@@#&~Nest~##&

Usage:

In[1]:= {+##/2,(1##)^.5}&@@#&~Nest~##&[{24, 6}, 5]

Out[1]= {13.4582, 13.4582}

Answer 18

R, 66 bytes

f=function(a,b,n){while(n){x=(a+b)/2;b=(a*b)^.5;n=n-1;a=x};c(a,b)}

Usage:

> f(24,6,0)
[1] 24  6
> f(24,6,1)
[1] 15 12
> f(24,6,2)
[1] 13.50000 13.41641
> f(24,6,3)
[1] 13.45820 13.45814
> f(24,6,4)
[1] 13.45817 13.45817
> f(100,50,10)
[1] 72.83955 72.83955
> f(1,1.41421356237,10)
[1] 1.19814 1.19814

Answer 19

Lua, 62 bytes

n,a,b=...for i=1,n do a,b=(a+b)/2,math.sqrt(a*b)end print(a,b)

Uses command line arguments from ... to assign to n, a and b, a nifty trick I learned about Lua recently.

Answer 20

Python 2, 62 61 62 bytes

def f(a,b,n):
 while n:a,b=(a+b)/2.,(a*b)**.5;n-=1
 print a,b

Answer 21

Haskell, 47 Bytes

f a b 0=(a,b)
f a b n=f((a+b)/2)(sqrt$a*b)(n-1)

Answer 22

Haskell, 40 bytes

(!!).iterate(\(a,b)->((a+b)/2,sqrt$a*b))

An anonymous function. Example usage:

>> let f=(!!).iterate(\(a,b)->((a+b)/2,sqrt$a*b)) in f (1.0,1.41421356237) 10
(1.198140234734168,1.1981402347341683)

The lambda function (\(a,b)->((a+b)/2,sqrt$a*b)) takes the arithmetic and geometric mean on a tuple. This is iterated starting with the first input (a tuple), and then (!!) indexes the second input to specify the number of iterations.

Answer 23

Perl, 60 bytes

perl -ape'F=($F[0]/2+$F[1]/2,sqrt$F[0]*$F[1])for 1..shift@F;$_="@F"'

N.B.: Per this meta post, I believe I've got the scoring correct. The actual code (between single quotes) is 58 characters, then I added +2 for a and p flags as that's the difference from the shortest invocation, perl -e'...'

Vague complaints

I have this nagging feeling I'm missing an obvious improvement. I know, "welcome to code golf", but I mean more than usual I believe there's an easy opportunity to shorten this.

Early on, I had messed around with using $\ as the second term with some success, but the above approach ended up being 2 bytes shorter, even with the extra ap flags required. Similarly, avoiding the explicit $_ assignment would be nice, but the loop makes that difficult.

The shift@F bugs me, too; if I don't do it that way, though (or use @F=(0,...,...) instead, which doesn't save any bytes), there's an off-by-one error with the @F assignment.

Example

echo 5 24 6 | perl -ape'F=($F[0]/2+$F[1]/2,sqrt$F[0]*$F[1])for 1..shift@F;$_="@F"'

Outputs

13.4581714817256 13.4581714817256

Answer 24

Julia, 49 bytes

(a,b,n)->(for i=1:n;a,b=(a+b)/2,√(a*b)end;(a,b))

Pretty direct iterative algorithm. Using the√ symbol and the multiple return saves a few bytes, but the for loop syntax costs a few.

Answer 25

Julia, 42 bytes

f(a,b,n)=n>0?f((a+b)/2,(a*b)^.5,n-1):(a,b)

This is a recursive function f that accepts three numbers and returns a tuple.

Ungolfed:

function f(a::Real, b::Real, n::Integer)
    if n > 0
        # Recurse on the arithmetic and geometric means, decrementing n
        return f((a + b) / 2, sqrt(a * b), n - 1)
    else
        # Return the pair
        return (a, b)
    end
end

Answer 26

LabVIEW, 21 LabVIEW Primitives

Primitives counted as per this meta post.

pretty staightforward not much to explain.

Answer 27

CJam, 16 bytes

{{_:+2/\:*mq]}*}

This is an anonymous function. The input is a list with the two values (as doubles), followed by the iteration count. Try it online with I/O code for testing.

I wouldn't normally have posted this because @PeterTaylor posted an equally long CJam answer before I saw the question. But since this is advertised as the start of a series, I wanted to keep my options open in case the series is interesting.

While the length is the same as Peter's answer, the code is not. I chose a different input format by taking the two values in a list, where Peter used separate values. So while there's not much to it with either input format, the code looks quite different.

{     Start loop over number of iterations.
  _     Copy the current pair of values.
  :+    Reduce pair with + operator.
  2/    Divide by 2.
  \     Swap second copy of pair to top.
  :*    Reduce pair with * operator.
  mq    Calculate square root.
  ]     Wrap the two new values in a list for next iteration.
}*    End iteration loop.

Answer 28

Perl 6,  53  47 bytes

{(($^a,$^b),->(\a,\b){((a+b)/2,sqrt(a*b))}...*)[$^n]} # 53 bytes

usage:

# give it a name
my &code = {(($^a,$^b),->(\a,\b){((a+b)/2,sqrt(a*b))}...*)[$^n]}

say code 100,50,10;          # (72.8395515523453 72.8395515523453)
say code 1,1.41421356237,10; # (1.19814023473417 1.19814023473417)

---

If I change the input from a,b,n to (a,b),n I can save a few bytes.

{($^l,->(\a,\b){((a+b)/2,sqrt(a*b))}...*)[$^n]} # 47 bytes

usage:

my &code = {($^l,->(\a,\b){((a+b)/2,sqrt(a*b))}...*)[$^n]}

say code (100,50),10;          # (72.8395515523453 72.8395515523453)
say code (1,1.41421356237),10; # (1.19814023473417 1.19814023473417)

say code (24,6),$_ for 0,1,2,5;
# (24 6)
# (15 12)
# (13.5 13.4164078649987)
# (13.4581714817256 13.4581714817256)

{
  (
    $^l,          # first 2 element tuple
    ->            # pointy block (lambda)
      (\a,\b)     # take a single tuple, and give its 2 elements each a name
    {
      (           # create a 2 element tuple
        (a+b)/2,  # arithmetic mean
        sqrt(a*b) # geometric mean
      )
    } ... *       # create a lazy infinite sequence of tuples
  )[ $^n ]        # take the nth "tuple" from the outer sequence
}

---

Really I would swap out the ... * with ... -> (\a,\b) { a =~= b }, then there would be no need for the $^n parameter.
( do not use == instead of =~=, or it may not stop )

my &code = {($^l,->(\a,\b){((a+b)/2,sqrt(a*b))}...->(\a,\b){a=~=b})[*-1]}

say code (24,6);           # (13.4581714817256 13.4581714817256)
say code (100,50);         # (72.8395515523453 72.8395515523453)
say code (1,1.41421356237) # (1.19814023473417 1.19814023473417)

Answer 29

Prolog, 80 bytes

Code:

p(A,B,0):-write([A,B]).
p(A,B,N):-X is(A+B)/2,Y is sqrt(A*B),M is N-1,p(X,Y,M).

Example:

p(100,50,10).
[72.83955155234534, 72.83955155234534]

Try it online here

Answer 30

Java, 103 96 84 bytes

String f(int n,double a,double b){return n>0?f(n-1,(a+b)/2,Math.sqrt(a*b)):a+","+b;}

Verify all testcases.

Old version (96 bytes):

String f(int n,double a,double b){for(;n>0;a=(a+b)/2,b=Math.sqrt((b-2*a)*b))n--;return a+","+b;}

Old version (103 bytes):

String f(int n,double a,double b){double t;for(;n>0;t=(a+b)/2,b=Math.sqrt(a*b),a=t)n--;return a+","+b;}