## Geothmetic Meandian

### April 13, 2021

In mathematics, the arithmetic geometric mean of two positive real numbers is computed by repeatedly taking half their sum and the square root of their product until the two numbers converge. For instance, the arithmetic geometric mean of 24 and 6 is 13.456171…, with the iterative steps computed as follows:

```0  24                            6
1  15                           12
2  13.5                         13.416407864998738175455042
3  13.458203932499369089227521  13.458139030990984877207090
4  13.458171481745176983217305  13.458171481706053858316334
5  13.458171481725615420766820  13.458171481725615420766806```

The arithmetic geometric mean was invented by Lagrange and studied by Gauss, and is used today to compute various transcendental functions because it converges so quickly.

In the world of Randall Munroe, the geothmetic meandian of any set of positive numbers is computed by iterating three sequences — the arithmetic mean, the geometric mean, and the median — until they converge. For instance, the geothmetic meandian of the set (1,1,2,3,5) is 2.089, computed as follows:

```1  2.4                1.9743504858348200 2
2  2.1247834952782734 2.1161924605448084 2
3  2.0803253186076938 2.0795368194795802 2.1161924605448084
4  2.0920181995440275 2.0919486049152223 2.0803253186076938
5  2.0880973743556477 2.0880901331209600 2.0919486049152223
6  2.0893787041306098 2.0893779142184865 2.0880973743556477
7  2.0889513309015815 2.0889512436159920 2.0893779142184865
8  2.0890934962453533 2.0890934865653277 2.0889513309015815
9  2.0890461045707540 2.0890461034958396 2.0890934865653277```

[ I hate the damned new editor at WordPress; I struggle with it every time I post an exercise. I could not figure out how to embed the image from XKCD in this blog post. You can see it here. ]

Your task is to write programs that compute the arithmetic geometric mean and geothmetic meandian. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

### 7 Responses to “Geothmetic Meandian”

1. Zack said

Fun little drill. I’m not sure about the usefulness of these centrality metrics, but it’s good practice though! Here is my take on this drill using Julia 1.6: https://pastebin.com/sj1X0G0H

2. pjbinformatimagocom said
```(defun arithmetic-geometric-mean (a b)
(values (/ (+ a b) 2)
(sqrt (* a b))))

(defun converge (epsilon fun a b)
(loop
:while (< epsilon (abs (- a b)))
:do (multiple-value-setq (a b) (funcall fun a b))
:finally (return (values a b))))

(converge 1d-15
(lambda (a b)
(format t ";; ~24,21F ~24,21F~%" a b)
(arithmetic-geometric-mean a b))
24d0 6d0)

;; 24.000000000000000000000  6.000000000000000000000
;; 15.000000000000000000000 12.000000000000000000000
;; 13.500000000000000000000 13.416407864998739000000
;; 13.458203932499370000000 13.458139030990985000000
;; 13.458171481745177000000 13.458171481706055000000
13.458171481725616D0
13.458171481725616D0

(defun geothmetic-meandian (set)
(sort (vector (/ (reduce (function +) set :initial-value 0.0d0)
(length set))
(expt (reduce (function *) set  :initial-value 1.0d0)
(/ (length set)))
(aref set (truncate (length set) 2)))
(function <)))

(defun converge-set (epsilon fun set)
;; set is sorted.
(loop
:while (< epsilon (abs (- (aref set 0) (aref set (1- (length set))))))
:do (setf set (funcall fun set))
:finally (return set)))

(converge-set
1d-3
(lambda (set)
(format t ";;~{ ~24,21F~}~%" (coerce set 'list))
(geothmetic-meandian set))
#(1 1 2 3 5))

;;  1.000000000000000000000  1.000000000000000000000  2.000000000000000000000  3.000000000000000000000  5.000000000000000000000
;;  1.974350485834820000000  2.000000000000000000000  2.400000000000000000000
;;  2.000000000000000000000  2.116192460544808400000  2.124783495278273000000
;;  2.079536819479580200000  2.080325318607693800000  2.116192460544808400000
;;  2.080325318607693800000  2.091948604915222300000  2.092018199544027500000
;;  2.088090133120960000000  2.088097374355647700000  2.091948604915222300000
;;  2.088097374355647700000  2.089377914218486500000  2.089378704130609800000
#(2.088951243615992D0 2.0889513309015815D0 2.0893779142184865D0)
```
3. pjbinformatimagocom said

You can run it at: https://ideone.com/8s3uh4

4. Daniel said

Here’s a solution in Python, utilizing the `statistics` module added to the standard library in Python 3.4.

```from statistics import geometric_mean, mean, median

EPSILON = 1e-10

def geometric_meandian(num1, num2):
while abs(num1 - num2) > EPSILON:
num1, num2 = mean((num1, num2)), geometric_mean((num1, num2))
return num1

def geothmetic_meandian(*nums):
while max(nums) - min(nums) > EPSILON:
nums = mean(nums), geometric_mean(nums), median(nums)
return nums

print(geometric_meandian(24, 6))
print(geothmetic_meandian(1, 1, 2, 3, 5))
```

Output:

```13.458171481745177
2.0890579497145714
```
5. Daniel said

Here’s another Python solution, which doesn’t rely on the `statistics` module.

```import math

EPSILON = 1e-10

def geometric_meandian(num1, num2):
while abs(num1 - num2) > EPSILON:
num1, num2 = num1 / 2 + num2 / 2, math.sqrt(num1 * num2)
return num1

def geothmetic_meandian(*nums):
nums = sorted(nums)
while nums[-1] - nums > EPSILON:
amean = sum(x / len(nums) for x in nums)
gmean = math.prod(nums) ** (1 / len(nums))
median = nums[len(nums) // 2]
if len(nums) & 1 == 0:
median = median / 2 + nums[len(nums) // 2 - 1] / 2
nums = sorted([amean, gmean, median])
return nums

print(geometric_meandian(24, 6))
print(geothmetic_meandian(1, 1, 2, 3, 5))
```

Output:

```13.458171481745177
2.0890579497145705
```
6. MSoni said

I have written the code of the following problem in Fortran using subroutine, so that it can be used for finding the solution of other inputted numbers also
<a href=”https://pastebin.com/ZQ7inPa4”>

7. Globules said

Here’s a Haskell version. I used this as an excuse to play around with the foldl library.

```import Control.Applicative (liftA2)
import qualified Control.Foldl as F
import Data.Bool (bool)
import Data.List (sort)
import Data.List.NonEmpty (NonEmpty, fromList, toList)
import Data.Maybe (fromJust)

{- ----------------------- Arithmetic-Geometric Mean ------------------------ -}

agm :: (Floating a, Ord a) => a -> (a, a) -> a
agm tol = fst . while ((> tol) . diff) step
where step (a, g) = ((a + g) / 2, sqrt (a * g))
diff (x, y) = abs (x - y)

{- --------------------------------- Median --------------------------------- -}

medianF :: (Fractional a, Ord a) => F.Fold a a
medianF = F.fold F.mean . middle <\$> sortF
where sortF = (,) <\$> F.length <*> (sort <\$> F.list)
middle (n, xs) = take (bool 1 2 (even n)) \$ drop ((n - 1) `div` 2) xs

{- ----------------------------- Geometric Mean ----------------------------- -}

gmeanF :: Floating a => F.Fold a a
gmeanF = nthRoot <\$> F.product <*> F.length
where nthRoot x n = x ** (1 / fromIntegral n)

{- -------------------------- Geothmetic Meandian --------------------------- -}

gmm :: (Floating a, Ord a) => a -> NonEmpty a -> a
gmm tol = head . gmms . toList
where gmms = while ((> tol) . maxDiff) (F.fold gmmF)
gmmF = sequenceA [F.mean, gmeanF, medianF]
maxDiff = fromJust . F.fold (liftA2 (-) <\$> F.maximum <*> F.minimum)

{- -------------------------------------------------------------------------- -}

while :: (a -> Bool) -> (a -> a) -> a -> a
while p f = head . dropWhile p . drop 1 . iterate f

test :: Show a => String -> (a -> Double) -> a -> IO ()
test name f arg = putStrLn \$ name ++ ": " ++ show arg ++ " => " ++ show (f arg)

main :: IO ()
main = do
let tol = 0.00001
xy = (24, 6)
ys = fromList [1, 1, 2, 3, 5]
test "AGM" (agm tol) xy
test "GMM" (gmm tol) ys
```
```\$ ./geomethic
AGM: (24.0,6.0) => 13.458171481745177
GMM: 1.0 :| [1.0,2.0,3.0,5.0] => 2.0890566336242085
```