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[0]
    
    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[0] > 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[0]
    
    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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: