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.
Programming Praxis 
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
(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)You can run it at: https://ideone.com/8s3uh4
Here’s a solution in Python, utilizing the
statisticsmodule 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:
Here’s another Python solution, which doesn’t rely on the
statisticsmodule.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:
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”>
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