Jump to content

Arithmetic-geometric mean

From Rosetta Code
Task
Arithmetic-geometric mean
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Arithmetic-geometric mean. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)


task

Write a function to compute the arithmetic-geometric mean of two numbers. [1] The arithmetic-geometric mean of two numbers can be (usefully) denoted as , and is equal to the limit of the sequence:

Since the limit of tends (rapidly) to zero with iterations, this is an efficient method.

Demonstrate the function by calculating:


Also see



360 Assembly

For maximum compatibility, this program uses only the basic instruction set. <lang 360asm>AGM CSECT 

        USING  AGM,R13

SAVEAREA B STM-SAVEAREA(R15) 

        DC     17F'0'
        DC     CL8'AGM'

STM STM R14,R12,12(R13) 

        ST     R13,4(R15)
        ST     R15,8(R13)
        LR     R13,R15
        ZAP    A,K                a=1
        ZAP    PWL8,K
        MP     PWL8,K
        DP     PWL8,=P'2'
        ZAP    PWL8,PWL8(7)
        BAL    R14,SQRT
        ZAP    G,PWL8             g=sqrt(1/2)

WHILE1 EQU * while a!=g 

        ZAP    PWL8,A
        SP     PWL8,G
        CP     PWL8,=P'0'         (a-g)!=0
        BE     EWHILE1
        ZAP    PWL8,A
        AP     PWL8,G
        DP     PWL8,=P'2'
        ZAP    AN,PWL8(7)         an=(a+g)/2
        ZAP    PWL8,A
        MP     PWL8,G
        BAL    R14,SQRT
        ZAP    G,PWL8             g=sqrt(a*g)
        ZAP    A,AN               a=an
        B      WHILE1

EWHILE1 EQU * 

        ZAP    PWL8,A
        UNPK   ZWL16,PWL8
        MVC    CWL16,ZWL16
        OI     CWL16+15,X'F0'
        MVI    CWL16,C'+'
        CP     PWL8,=P'0'
        BNM    *+8
        MVI    CWL16,C'-'
        MVC    CWL80+0(15),CWL16
        MVC    CWL80+9(1),=C'.'   /k  (15-6=9)
        XPRNT  CWL80,80           display a
        L      R13,4(0,R13)
        LM     R14,R12,12(R13)
        XR     R15,R15
        BR     R14
        DS     0F

K DC PL8'1000000' 10^6 A DS PL8 G DS PL8 AN DS PL8

  • ****** SQRT *******************

SQRT CNOP 0,4 function sqrt(x) 

        ZAP    X,PWL8
        ZAP    X0,=P'0'           x0=0
        ZAP    X1,=P'1'           x1=1

WHILE2 EQU * while x0!=x1 

        ZAP    PWL8,X0
        SP     PWL8,X1
        CP     PWL8,=P'0'         (x0-x1)!=0
        BE     EWHILE2
        ZAP    X0,X1              x0=x1
        ZAP    PWL16,X
        DP     PWL16,X1
        ZAP    XW,PWL16(8)        xw=x/x1
        ZAP    PWL8,X1
        AP     PWL8,XW
        DP     PWL8,=P'2'
        ZAP    PWL8,PWL8(7)
        ZAP    X2,PWL8            x2=(x1+xw)/2
        ZAP    X1,X2              x1=x2
        B      WHILE2

EWHILE2 EQU * 

        ZAP    PWL8,X1            return x1
        BR     R14
        DS     0F

X DS PL8 X0 DS PL8 X1 DS PL8 X2 DS PL8 XW DS PL8

  • end SQRT

PWL8 DC PL8'0' PWL16 DC PL16'0' CWL80 DC CL80' ' CWL16 DS CL16 ZWL16 DS ZL16 

        LTORG  
        YREGS  
        END    AGM</lang>
Output:
+00000000.84721

8th

<lang 8th>: epsilon 1.0e-12 ;

with: n

iter \ n1 n2 -- n1 n2
   2dup * sqrt >r + 2 / r> ;
agn \ n1 n2 -- n
   repeat  iter  2dup epsilon ~= not while!  drop ;

"agn(1, 1/sqrt(2)) = " . 1 1 2 sqrt / agn "%.10f" s:strfmt . cr

with

bye </lang>

Output:
agn(1, 1/sqrt(2)) = 0.8472130848

Ada

<lang Ada>with Ada.Text_IO, Ada.Numerics.Generic_Elementary_Functions;

procedure Arith_Geom_Mean is 

  type Num is digits 18; -- the largest value gnat/gcc allows
  package N_IO is new Ada.Text_IO.Float_IO(Num);
  package Math is new Ada.Numerics.Generic_Elementary_Functions(Num);

  function AGM(A, G: Num) return Num is
     Old_G: Num;
     New_G: Num := G;
     New_A: Num := A;
  begin
     loop
        Old_G := New_G;
        New_G := Math.Sqrt(New_A*New_G);
        New_A := (Old_G + New_A) * 0.5;
        exit when (New_A - New_G) <= Num'Epsilon;
        -- Num'Epsilon denotes the relative error when performing arithmetic over Num
     end loop;
     return New_G;
  end AGM;

begin 

  N_IO.Put(AGM(1.0, 1.0/Math.Sqrt(2.0)), Fore => 1, Aft => 17, Exp => 0);

end Arith_Geom_Mean;</lang>

Output:

0.84721308479397909

ALGOL 68

Algol 68 Genie gives IEEE double precision for REAL quantities, 28 decimal digits for LONG REALs and, by default, 63 decimal digits for LONG LONG REAL though this can be made arbitrarily greater with a pragmat.

Printing out the difference between the means at each iteration nicely demonstrates the quadratic convergence. <lang algol68> BEGIN 

  PROC agm = (LONG REAL x, y) LONG REAL :
  BEGIN
     IF x < LONG 0.0 OR y < LONG 0.0 THEN -LONG 1.0
     ELIF x + y = LONG 0.0 THEN LONG 0.0		CO Edge cases CO
     ELSE

LONG REAL a := x, g := y; LONG REAL epsilon := a + g; LONG REAL next a := (a + g) / LONG 2.0, next g := long sqrt (a * g); LONG REAL next epsilon := ABS (a - g); WHILE next epsilon < epsilon DO print ((epsilon, " ", next epsilon, newline)); epsilon := next epsilon; a := next a; g := next g; next a := (a + g) / LONG 2.0; next g := long sqrt (a * g); next epsilon := ABS (a - g) OD; a 

     FI
  END;
  printf (($l(-35,33)l$, agm (LONG 1.0, LONG 1.0 / long sqrt (LONG 2.0))))

END </lang>

Output:

+1.707106781186547524400844362e  +0   +2.928932188134524755991556379e  -1
+2.928932188134524755991556379e  -1   +1.265697533955921916929670477e  -2
+1.265697533955921916929670477e  -2   +2.363617660269221214237489508e  -5
+2.363617660269221214237489508e  -5   +8.242743980540458935740117000e -11
+8.242743980540458935740117000e -11   +1.002445937606580000000000000e -21
+1.002445937606580000000000000e -21   +4.595001000000000000000000000e -29
+4.595001000000000000000000000e -29   +4.595000000000000000000000000e -29

0.847213084793979086606499123550000

APL

<lang APL> agd←{(⍺-⍵)<10*¯8:⍺⋄((⍺+⍵)÷2)∇(⍺×⍵)*÷2} 1 agd ÷2*÷2 </lang>

Output: 

0.8472130848

AppleScript

By functional composition:

<lang AppleScript>-- ARITHMETIC GEOMETRIC MEAN -------------------------------------------------

property tolerance : 1.0E-5

-- agm :: Num a => a -> a -> a on agm(a, g) 

   script withinTolerance
       on |λ|(m)
           tell m to ((its an) - (its gn)) < tolerance
       end |λ|
   end script
   
   script nextRefinement
       on |λ|(m)
           tell m
               set {an, gn} to {its an, its gn}
               {an:(an + gn) / 2, gn:(an * gn) ^ 0.5}
           end tell
       end |λ|
   end script
   
   an of |until|(withinTolerance, ¬
       nextRefinement, {an:(a + g) / 2, gn:(a * g) ^ 0.5})

end agm

-- TEST ---------------------------------------------------------------------- on run 

   agm(1, 1 / (2 ^ 0.5))
   
   --> 0.847213084835

end run

-- GENERIC FUNCTIONS ---------------------------------------------------------

-- until :: (a -> Bool) -> (a -> a) -> a -> a on |until|(p, f, x) 

   set mp to mReturn(p)
   set v to x
   tell mReturn(f)
       repeat until mp's |λ|(v)
           set v to |λ|(v)
       end repeat
   end tell
   return v

end |until|

-- Lift 2nd class handler function into 1st class script wrapper -- mReturn :: Handler -> Script on mReturn(f) 

   if class of f is script then
       f
   else
       script
           property |λ| : f
       end script
   end if

end mReturn</lang>

Output:
0.847213084835

AutoHotkey

<lang AHK>agm(a, g, tolerance=1.0e-15){ While abs(a-g) > tolerance { an := .5 * (a + g) g  := sqrt(a*g) a  := an } return a } SetFormat, FloatFast, 0.15 MsgBox % agm(1, 1/sqrt(2))</lang> Output: 

0.847213084793979


AWK

<lang AWK>#!/usr/bin/awk -f BEGIN { printf "%.16g\n", agm(1.0,sqrt(0.5)) } function agm(a,g) { while (1) { a0=a a=(a0+g)/2 g=sqrt(a0*g) if (abs(a0-a) < abs(a)*1e-15) break } return a } function abs(x) { return (x<0 ? -x : x) } </lang> Output 

0.8472130847939792

BASIC

Commodore BASIC

<lang commodorebasic>10 A = 1 20 G = 1/SQR(2) 30 GOSUB 100 40 PRINT A 50 END 100 TA = A 110 A = (A+G)/2 120 G = SQR(TA*G) 130 IF A<TA THEN 100 140 RETURN</lang>


BBC BASIC

Works with: BBC BASIC for Windows

<lang bbcbasic> *FLOAT 64 

     @% = &1010
     PRINT FNagm(1, 1/SQR(2))
     END
     
     DEF FNagm(a,g)
     LOCAL ta
     REPEAT
       ta = a
       a = (a+g)/2
       g = SQR(ta*g)
     UNTIL a = ta
     = a

</lang> Produces this output: 

0.8472130847939792

bc

<lang bc>/* Calculate the arithmethic-geometric mean of two positive 

* numbers x and y.
* Result will have d digits after the decimal point.
*/

define m(x, y, d) { 

   auto a, g, o

   o = scale
   scale = d
   d = 1 / 10 ^ d

   a = (x + y) / 2
   g = sqrt(x * y)
   while ((a - g) > d) {
       x = (a + g) / 2
       g = sqrt(a * g)
       a = x
   }

   scale = o
   return(a)

}

scale = 20 m(1, 1 / sqrt(2), 20)</lang>

Output:
.84721308479397908659

C

Basic

<lang c>#include<math.h>

  1. include<stdio.h>
  2. include<stdlib.h>

double agm( double a, double g ) { 

  /* arithmetic-geometric mean */
  double iota = 1.0E-16;
  double a1, g1;
  
  if( a*g < 0.0 ) {
     printf( "arithmetic-geometric mean undefined when x*y<0\n" );
     exit(1);
  }

  while( fabs(a-g)>iota ) {
     a1 = (a + g) / 2.0;
     g1 = sqrt(a * g);

     a = a1;
     g = g1;
  }
  
  return a;

}

int main( void ) { 

  double x, y;
  printf( "Enter two numbers: " );
  scanf( "%lf%lf", &x, &y );
  printf( "The arithmetic-geometric mean is %lf\n", agm(x, y) );
  return 0;

} </lang>

Original output: 

Enter two numbers: 1.0 2.0
The arithmetic-geometric mean is 1.456791

Task output, the second input (0.707) is 1/sqrt(2) correct to 3 decimal places: 

Enter two numbers: 1 0.707
The arithmetic-geometric mean is 0.847155

GMP

<lang cpp>/*Arithmetic Geometric Mean of 1 and 1/sqrt(2) 

 Nigel_Galloway
 February 7th., 2012.
  • /
  1. include "gmp.h"

void agm (const mpf_t in1, const mpf_t in2, mpf_t out1, mpf_t out2) { mpf_add (out1, in1, in2); mpf_div_ui (out1, out1, 2); mpf_mul (out2, in1, in2); mpf_sqrt (out2, out2); }

int main (void) { mpf_set_default_prec (65568); mpf_t x0, y0, resA, resB;

mpf_init_set_ui (y0, 1); mpf_init_set_d (x0, 0.5); mpf_sqrt (x0, x0); mpf_init (resA); mpf_init (resB);

for(int i=0; i<7; i++){ agm(x0, y0, resA, resB); agm(resA, resB, x0, y0); } gmp_printf ("%.20000Ff\n", x0); gmp_printf ("%.20000Ff\n\n", y0);

return 0; }</lang>

The first couple of iterations produces: 

0.853
0.840

Then 7 iterations produces: 

0.84721308479397908660649912348219163648144591032694218506057937265973400483413475972320029399461122994212228562523341096309796266583087105969971363598338425117632681428906038970676860161665004828118872189771330941176746201994439296290216728919449950723167789734686394760667105798055785217314034939830420042211921603983955359509819364129371634064602959996797059943435160203184264875695024217486385540598195458160174241788785419275880416271901208558768564832683414043121840080403580920455949431387781512092652225457439712428682076634095473367459962179266553534862568611854330862628728728756301083556319357066871478563908898211510883635214769697961262183294322841786811376844517001814602191369402702094599668351359632788080427434548174458736322002515395293626580661419836561649162625960743472370661690235308001737531284785255843063190745427493415268579065526940600314759102033274671968612479632551055464890282085529743965124994009662552866067580448735389218570140116771697653501408495247684899325732133702898466893919466586187375296638756226604591477704420468108925658440838032040910619003153706734119594101007474331059905505820524326009951692792417478216976781061683697714110739273343921550143022007087367365962272149258776192851052380367026890463909621907663644235538085902945234065190013342345105838341712180514255003923701111325411144612628906254133550526643653595824552156293397518251470650134641047056979355681306606329373345038710977097294875917179015817320281578288487149931340815493342367797044712785937618595085146677364554679201615934223997142984070788882279032656751596528435817795727284808356489963504404140734226110183383546975962663330422084999852300742703930277243474979717973264552546543019831694968461098690743905068013766119252919770938441299707015889493166661161994592265011311183966352502530561646431587208454522988775475177272747656721648982918239238895207207642839710884705960356921992921831901548141280766592698294464457149239666329973075813904957622438962423175209507319018424462442370986427281149511180822826053862484617675180140983127497257651983756492356902800216174905531427208153439540595563576371127281657059737337442970039056040156388663072225700389230159112376960121580081779077863351240862431073571583765926504546652787337874444834406310244757039681255453982266430353416413035613801634165575265589752944521166873451220191227466733191571240763753821106968141076926390074833175743396752319660330864973571383874196098983832202882694882191302819366949954422240697276168621369511657838885012199096160655454611543253148164249332694797004159491476323112920593516518997943350045976288217292625918089405508431466393782548335139550190653370872062064024077056075848796499843651592728264534428636615419142585777106756185017278033287175195189305031805505245426022335522900771418128798654351187918006356279593624768267786412249460338126082628254098895312527677534656243279214511229555516031818433133692961723041783855157125567404983416665926969580008953724573057694542275372160209687191470398878466367243262706191127071716590824640041679941120405657103640830002419294398553073994656539677810492701055410359513339432199925066676202078394695553760551796401009749218856311301017813888578793813172095948062539201300983650287917695827985905279947721941797997024943062158419468885328115497721579960194409623477686144085075739284298823759396823223670580334134774623112897625859324376631778974911077261909704489522204509630725515590093824904021364807792034767215048568446022554409992826163174312642285787628983380650722023010371753149263504631060188573772567006618381290580638954508127031311371043716135833488065833955431217901348398833216413057635244712511539472066670330101348716516324113828817639839629526121141263219795965098656786755250760760424095907517523021946104532564333249614901253533329223723868948127885020135966305376055849358928391630469403887854960027471487197801457659579049585802260066099524967364324966833461760106608156706975142381866503610838852209761655002516073114992161294775790199729248689638220603808760276281672370166819106633585775154650381334236722347642026558565588464160102105404898556187114735884976378406486426798186504486319077470382286711435151123003607086574298864771466747337501143458188527970060562117246921748471806948662511994728934442703783046207073549380528727206215606307188286858056452111069670802856990698257691772209986719599685077906814434949328049768115436804632599386930762350709995182951295811212357072453833548261907523951582730982481805496658979091688679840717077937059590457758409104734131096041941113577566207273377978332037973011376726585357477102797814097213096121423938547374627696150413079528373728820506587191522597650840277969917611753930067254924912298450823629755687227110658494355338504945326387364898046066559799543601695030927900924500578564772358761988489860344121953407953690029964119745490607416009788595376607229051607724285900709011566391383642990412208267696297978676490323564999819907659974398705486487690910249119270999682756970113687622440464029603837000662127345776647097113263746568115029858630322603373834213584239378961146171920830719539156437820936414967803341524645073966831731983633627433925553117120194541468448808956224178980318943412312840278583782890096242095413450021010727363232852725762096468519944682405506293917420533017064619172151788442967053143355037723107097160802851453141441061050231173108777799332489320877272298978213301208340743056049981599632026877933071569403024391561189267675172495117665262485470960419914731136579206973309960888972867897807355875785006235751571237716530420636310027031292966940254219678771688466557275808983064676620070146795856930822206209053308277822265031125202787335125191599188939002843192181666865484348796219722117639049598957936073309436974576289432003841175529415947547471839363811441256103510234595810807685589856570074453089094286692511901017181228266893492695282610525185567360458777022881478214469685009183472197414205461280723479500598117663645261501907885454711938035571459307446356562607527875188243864095069646498151311705914579906193765608586501756168645019240983272357243336888130800221863687002096411197243036035586497937733149167495931511886735350255059823030470602847404584566768496209345063963029094416325164086928898145072478777276733780338289295049783843429437665667372975874305751410364174768616396241989419047309961002284280794449200269048452541391882460015590891319432556103657693623641617846466931414561099840383122655041152514944453800420904287181824684316246105526376775209701040639446878373750174360897516934868876512834536775527865470902315420294538730761411966497675219198089021057726334723979589687229233577690412444586822978062098870898160181795214549203709562528507330232550600966113294791484434166874298726542040835520564564044211741240650419323628312966431263307687154504449507335544182007936697013312446388243600624398167124093468063221697717015635904176098412619778010525869566346541447025111353828410102785795430618023572755009305139556377710439227995971141182782033581183989523387201196266668287812153433311933530198006525119241035943150724272515897742269014313251497752206211486532095282917841726788527918259501894283066454533808294385484913906600901526463156669408130516898577384457161101347735284395586639180314771289972489772326950830959208603163908601794221468048925371471356694906475975663504050761059303001534536134468346141362848404730639095800648624822113995399621221079927740532030597569871315014292389418219892184458614968453063460782870588642625603497671133853907530473607475205697255326635179640594881381276485191302328261295517207475944988639251110497859774101046472588317449694892733322810684089494759787067690122169518696581944061366943103234116196131605543816087283055435048190711597527426659173636930019809887976272186626285433119060860342806191518452978237036398984494144178890086027822209983902274728379674114295789243465456404028551674783725388313861547805080352368935833328873558797948868049809714068689367194167115043074025751022690817073859285358373909764249759224210618323725170214283209867537445071332189636669085656349633060774556830118371494002584049977661135255328476656188705929782127298997295927947818204287198071022786461838070064010831389756771127541362211274445345355849597692525757583129990395369598932499513241067842656115567436600887374842740382348117849110021235371080153344077081752815794229285487316898639800718962686849857790619425820001731784737979758156092690872878502700244147412819535788739647458594598995355434128016535530490585287946743982206062303866888527005052189049277821975141155954355491253261150874322804356095631761163218117941648842069284743156991336777879569137055927049598939111007862241124499317195398903082153071269718073528142944373740581805897842871015663258737266000122961804037804290931751604739799312368824663145245907925120889169747654302457053206386704684110540342014376644422132127507998462991570101471065529461467463922495745306196822034254448162475459772696534302506868242052880996924489236521714038177492829359173154812849196214333040809043068672336820607162912893985174062559042822475581595091023242061608163635114409532679679744662146581218973837257052018318006785051812332707432360517602365653046059197282467620464979507571243323062106152366172293244682862511105778328547123718579064823024291991297534773406188123932244051237932292486982393020946057994685022093564580188647372057989508199682850879081206451754647928466570299934961463545338169898790120739595342994580518846829188356311361388796313161734422075062182129450475034337306401403566141064033208676214431839284389699942682868360825355912427514883833922646682229633236574889815991049023745712780770628532368956900284697429547742484223355238590492992254533182706939660886035184911668751085520062653409664126112200692905563690527440648936400870151716629293565299214744207938737106473991364534021859315182015761100594055566001663181909163482128186430684182569911943162667158985886736504889805808329721451958115258329743580644326982892093642849596169753399275023838326958011096089547864572561097853782973070749181687447357311890498494907816322101271109193983576388927531317499783213682809328943493309300878688841270920763590076480651183013174408131381707764785620869834568499576963332415566990859371495284373037821741667810126247377548449594082775980428578137754484461929295371533597418713555566780286064849179748275590223773761897037703324897743492353765235571390764314889671441330995396798710462847477217721858658519859712821657391485744943283203084641639560963010473704739884503079369569286834641137642263085686956881520537491962945628810859870159107649550192726673782765172374500136624210511467091848989522697276562069762630550949389320992163775294153350600271094300189773392218453903373510079427646652325090453779404782123556204886389696402910291826730243688880139827500496556889555403627397541183592770090942918399583962985359521234655737077516804320238724010087862923625584849202212960559482323176352142071176504276997478012902491509148733472049812083534865212462335388584717004701205923945825415223129676013072682802320446336442341000264743415683991238810480498194912009402448957203018812206409969973408437360958124499459132317933593338191973602488533756410304356437323020013283599906152983949167106879976939266990335220640837295869943043576709171697966984423326568307325500003213129027067191063424283113900494781793073045562199439120722094954719165471096054049199441860517249814718129940631192901737381011766173569764956366756202788955920995046861634403052506586817358402694287366334311678329038374756580509907839853849260647212465651306604876736085857902183866432416271982103787727963377367426929456639854705293777458546922070020463303573435055175370140503103555265780827298970492305475455890092754109445040141571253576828010749151746279285337830995706319528768382378063681778416611863347477894201661901861433888045148841743616814548103623210376432745956533646293972952940499526616911816577400181161464976544075891509125575991008552731077337032136035056194073504052234145332243066047436002572125901272025171469526054624392158151517326614548122436198603573869224654036885597877500832683869306742537593493769726913825327805701356834418623150103189551287054940385947609492785905200098814477158397147139718137205549603311916422391953132302138759927174019046224139259148006201715618158893529451219781937047457085386954279002330804105880072509475123189307968446372241711705946061976147519773238961013155564063723093102794769739382294763468939337559468936650940499102526121635380720056442410264711646398004909985355702820593960545544792555586249187092321801304541029363328936193265963508514136372072931427677632678178400667800895586548777826308228184465081585096256950206977978896641405511014211855334440159488802847016579044649263092161202380685664726316113269955335854143205474428967281732917140106437305939602224827339697208658091942888039633443448764675833855973513333306284397863570621963822177055006726076075702023055483284393359373696240854049573444151418891438122060768323290633843326859359282266483616228768156709313037896783277414878452878382324740383408934494278060455890181836731336022711672853044271945073157409136000663560891812190403050193190281639721357906960252119295624559528358504426277879932144682210413256122712903024696103748551345991066626060821435461264637908469523386805592378228286103613864160137539204268883711926027420874745077827301808826482979914892334346533639303279918164769955294688929040603354702651883178258213919150731170223368395649453356304141924428385039542090733375111170537908197680613788461570042923922647881382284866725434155806944211935068360004884655615990833391847242631836989281306956549491531650103132163612240182987115172224015233681014762461698964172597488387271895987656023503248287097414687934153787088145731903279204532192316858527351083720559424566015456479446754495668591429979882331798190595741253686810321947980826038762410448487302089050658719342641740920079366698836014623097627598441130715257589162880105817093530725888876543862532018486249319236385682165626031104345283130307049722913348730332409337369563479748898249300174158056591821232883438581012501715373053984620434324557214820885475234947304677614292829153914858526885054230744505481926191669759750315034472082118453139076834860069087727520772464857065976367409361731434369903994989083757102465456508149620159888052044833794917070408483039094175124262758698686686442934982424196674036270760323992014071830712707598371320007124471595236427821624884729339137136340461389740888941783993200900515436084216188913289577403543844561076450160104627095790986524953420147660163304582935376534545234386674137987312550170295545828095478975424973671090385982646068956222412573032081408906070252061404578152823685045057657100438042285920327207291902221346518359302559429408753069947011011534164767856235435750239937364145328957734998761675022409197941218931880590179774443294036240385510824919547518411770141508205549991488032865000650690301650284556165335148907119741941723100296632479366408253645421048976404451080811239063681885949086604183400256315626612115063653092972195806871776320514613555813095008145638261124165214871635936435536462688727462766803686306800882312499705727064962653352854242737234497574827760613008180634196390830978822494789229495258916657826100444244401103267485396201200233971298346242423632837110742673099021260291100381090507518405232662739050319348560154855106326243187789708788951981680730963542230960055362677359050994734087443710248167279700094945897076301853449526801067309842468288488837600166958871373559692445552385363961787881342093093764848484068429404997314946635784558266882458253566353932897293167000662381283685196706276978897699290095978380695574407690809500695946595783253660660602130005250129981452150996293071107006157960047599188298274727518774924726747707554136792657750601495283368598380853534208742156827588012599928559034100979630199437410013949755918229188467057410106349315945279547420320572953565968695868630973284883811742438270584417356596674853152028861911921252863987395609281275132232141197542293430923755693396146727405175695293766990610523654483440786104255766945418734863793560708612404736883567734371401263501208237651763905620506040768947294002931620797603428968468976398678305539415152307137255605029146711751234519321319625717919409117289511239481135988605880624240378357519964870883301506792101754290605314188369786110278968306896668518684104701823647807006155298831498831116019499658150386743904671052471759937267092033810519847770061227523026980385376199177319071331058167790086514801724404464037647206737845833953828893809029412739879104752542584865616980485432967822810404539976611651232907291616199926287510865193417311165133056591829817625847694287084548190293442221860279774055192912661889487080105159228601492383934908897821669651094997616731795835221057913587243550297821114252805843809597704721778938273829164718826714378658214613260112635165542805164184221882641418906866191864927517189847350374966026860336719613049159226094421467730920744767947119178202099132268721849475483780038487261488727428812655791747946341514445451055994645676144782933879680154128864180982848855259596173991776576352670819899854089307445641992969024592754051436475256486619329599030683238667575184797410153429114165087535728924796842802484402202118983902434301907465924705639919100242258143990683914578574580953440968261584897316158220398376910051716543905900933268275864197534394837719059730794650292103636419726159238721878760956871976819344819558525670241414336715908896942047817989365563517751015910050265859472794486423173118927271535250460340818962273831146005468524063988554718596840882777221622505863684193799641126463210706398187737943696502521044386223206715172284114754334828030417076754385554475843212718463962813919258849725090510409441344504298453460718488756542407096901385926116455196765637084297106764946357662012853819267912041109778058573520627375104669435915920749043789661298087162743223850390320074778542110638995449541859976414281163951972397080789860487582641265448251499232272861765713896973345378359636039627090380026689213243891590093752250336511719377706572262953412570689809077931988799970767832633036706673426579253958499505823639986104928784799761858913840247447907423559817960132549606526849887335183972871912518993883243416026083561644966709023900422732162219315679399440012151599100543810845200811331032075534924844873692683144444666107802758917774683693445850459499632371560438002582276189086030745508199318928997032855495073302401217663495153158278308977864322545562217443057528251437080871843144708110045101086121226999313969693610665236087211263590123448282622844271912819731872697619747403980717783781881605198018622572329702247624947679129326840201880617952362291746013985766042335790944077230173530153379744356437385842482505380615471930752244293091172074476771495221419193909742017160269705578258369237072978115455525707880049556669154779018307195916635166870579843369516111891537519123967141163781970007849531153863267663692691720169784090403969698048618284364177768040884492084399010959512057513408610603753534081557370871883138983376563225336509460103086861119012415417949006598353669263835150584020260982595703854291458650256921579873098070645970823263771382355857377042256281442627934977694293588040208827420282637864436159358179308178583062657122634794521740652164107980293335739611374043019282943678846268324324490788126847872819886762029310625102649485865494639647891543662406355703466884777848152714124704306460406156142773201070035758550339952793775297161566283811185180855234141875772560252179951036627714775522910368395398000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

The limit (19,740) is imposed by the accuracy (65568). Using 6 iterations would produce a less accurate result. At 7 iterations increasing the 65568 would mean we already have 38,000 or so digits accurate.

C++

<lang c++>

  1. include<bits/stdc++.h>

using namespace std;

  1. define _cin ios_base::sync_with_stdio(0); cin.tie(0);
  2. define rep(a, b) for(ll i =a;i<=b;++i)

double agm(double a, double g) //ARITHMETIC GEOMETRIC MEAN { double epsilon = 1.0E-16,a1,g1; if(a*g<0.0) { cout<<"Couldn't find arithmetic-geometric mean of these numbers\n"; exit(1); } while(fabs(a-g)>epsilon) { a1 = (a+g)/2.0; g1 = sqrt(a*g); a = a1; g = g1; } return a; }

int main() { _cin; //fast input-output double x, y; cout<<"Enter X and Y: "; //Enter two numbers cin>>x>>y; cout<<"\nThe Arithmetic-Geometric Mean of "<<x<<" and "<<y<<" is "<<agm(x, y); return 0; } </lang> 


Enter X and Y: 1.0 2.0 
The Arithmetic-Geometric Mean of 1.0 and 2.0 is 1.45679103104690677028543177584651857614517211914062


C#

<lang csharp>namespace RosettaCode.ArithmeticGeometricMean { 

   using System;
   using System.Collections.Generic;
   using System.Globalization;

   internal static class Program
   {
       private static double ArithmeticGeometricMean(double number,
                                                     double otherNumber,
                                                     IEqualityComparer<double>
                                                         comparer)
       {
           return comparer.Equals(number, otherNumber)
                      ? number
                      : ArithmeticGeometricMean(
                          ArithmeticMean(number, otherNumber),
                          GeometricMean(number, otherNumber), comparer);
       }

       private static double ArithmeticMean(double number, double otherNumber)
       {
           return 0.5 * (number + otherNumber);
       }

       private static double GeometricMean(double number, double otherNumber)
       {
           return Math.Sqrt(number * otherNumber);
       }

       private static void Main()
       {
           Console.WriteLine(
               ArithmeticGeometricMean(1, 0.5 * Math.Sqrt(2),
                                       new RelativeDifferenceComparer(1e-5)).
                   ToString(CultureInfo.InvariantCulture));
       }

       private class RelativeDifferenceComparer : IEqualityComparer<double>
       {
           private readonly double _maximumRelativeDifference;

           internal RelativeDifferenceComparer(double maximumRelativeDifference)
           {
               _maximumRelativeDifference = maximumRelativeDifference;
           }

           public bool Equals(double number, double otherNumber)
           {
               return RelativeDifference(number, otherNumber) <=
                      _maximumRelativeDifference;
           }

           public int GetHashCode(double number)
           {
               return number.GetHashCode();
           }

           private static double RelativeDifference(double number,
                                                    double otherNumber)
           {
               return AbsoluteDifference(number, otherNumber) /
                      Norm(number, otherNumber);
           }

           private static double AbsoluteDifference(double number,
                                                    double otherNumber)
           {
               return Math.Abs(number - otherNumber);
           }

           private static double Norm(double number, double otherNumber)
           {
               return 0.5 * (Math.Abs(number) + Math.Abs(otherNumber));
           }
       }
   }

}</lang> Output: 

0.847213084835193

Clojure

<lang lisp>(ns agmcompute 

 (:gen-class))
Java Arbitray Precision Library

(import '(org.apfloat Apfloat ApfloatMath))

(def precision 70) (def one (Apfloat. 1M precision)) (def two (Apfloat. 2M precision)) (def half (Apfloat. 0.5M precision)) (def isqrt2 (.divide one (ApfloatMath/pow two half))) (def TOLERANCE (Apfloat. 0.000000M precision))

(defn agm [a g] 

 " Simple AGM Loop calculation "
      (let [THRESH 1e-65                 ; done when error less than threshold or we exceed max loops
            MAX-LOOPS 1000000]
       (loop [[an gn] [a g], cnt 0]
           (if (or (< (ApfloatMath/abs (.subtract an gn)) THRESH)
                   (> cnt MAX-LOOPS))
             an
             (recur [(.multiply (.add an gn) half) (ApfloatMath/pow (.multiply an gn) half)]
                    (inc cnt))))))

(println (agm one isqrt2)) </lang>

Output:
8.47213084793979086606499123482191636481445910326942185060579372659734e-1

COBOL

<lang cobol>IDENTIFICATION DIVISION. PROGRAM-ID. ARITHMETIC-GEOMETRIC-MEAN-PROG. DATA DIVISION. WORKING-STORAGE SECTION. 01 AGM-VARS. 

   05 A       PIC 9V9(16).
   05 A-ZERO  PIC 9V9(16).
   05 G       PIC 9V9(16).
   05 DIFF    PIC 9V9(16) VALUE 1.
  • Initialize DIFF with a non-zero value, otherwise AGM-PARAGRAPH
  • is never performed at all.

PROCEDURE DIVISION. TEST-PARAGRAPH. 

   MOVE    1 TO A.
   COMPUTE G = 1 / FUNCTION SQRT(2).
  • The program will run with the test values. If you would rather
  • calculate the AGM of numbers input at the console, comment out
  • TEST-PARAGRAPH and un-comment-out INPUT-A-AND-G-PARAGRAPH.
  • INPUT-A-AND-G-PARAGRAPH.
  • DISPLAY 'Enter two numbers.'
  • ACCEPT A.
  • ACCEPT G.

CONTROL-PARAGRAPH. 

   PERFORM AGM-PARAGRAPH UNTIL DIFF IS LESS THAN 0.000000000000001.
   DISPLAY A.
   STOP RUN.

AGM-PARAGRAPH. 

   MOVE     A TO A-ZERO.
   COMPUTE  A = (A-ZERO + G) / 2.
   MULTIPLY A-ZERO BY G GIVING G.
   COMPUTE  G = FUNCTION SQRT(G).
   SUBTRACT A FROM G GIVING DIFF.
   COMPUTE  DIFF = FUNCTION ABS(DIFF).</lang>
Output:
0.8472130847939792

Common Lisp

<lang lisp>(defun agm (a0 g0 &optional (tolerance 1d-8)) 

 (loop for a = a0 then (* (+ a g) 5d-1)
    and g = g0 then (sqrt (* a g))
    until (< (abs (- a g)) tolerance)
    finally (return a)))

</lang>

Output:
CL-USER> (agm 1d0 (/ 1d0 (sqrt 2d0)))
0.8472130848351929d0
CL-USER> (agm 1d0 (/ 1d0 (sqrt 2d0)) 1d-10)
0.8472130848351929d0
CL-USER> (agm 1d0 (/ 1d0 (sqrt 2d0)) 1d-12)
0.8472130847939792d0

D

<lang d>import std.stdio, std.math, std.meta, std.typecons;

real agm(real a, real g, in int bitPrecision=60) pure nothrow @nogc @safe { 

   do {
       //{a, g} = {(a + g) / 2.0, sqrt(a * g)};
       AliasSeq!(a, g) = tuple((a + g) / 2.0, sqrt(a * g));
   } while (feqrel(a, g) < bitPrecision);
   return a;

}

void main() @safe { 

   writefln("%0.19f", agm(1, 1 / sqrt(2.0)));

}</lang>

Output:
0.8472130847939790866

All the digits shown are exact.

EchoLisp

We use the (~= a b) operator which tests for |a - b| < ε = (math-precision). <lang scheme> (lib 'math)

(define (agm a g) 

   (if (~= a g) a 
      (agm (// (+ a g ) 2) (sqrt (* a g)))))

(math-precision) 

   → 0.000001 ;; default

(agm 1 (/ 1 (sqrt 2))) 

   → 0.8472130848351929

(math-precision 1.e-15) 

   → 1e-15

(agm 1 (/ 1 (sqrt 2))) 

   → 0.8472130847939792

</lang>


Elixir

<lang Elixir>defmodule ArithhGeom do 

 def mean(a,g,tol) when abs(a-g) <= tol, do: a
 def mean(a,g,tol) do
   mean((a+g)/2,:math.pow(a*g, 0.5),tol)
 end

end

IO.puts ArithhGeom.mean(1,1/:math.sqrt(2),0.0000000001)</lang>

Output:
0.8472130848351929

Erlang

<lang Erlang>%% Arithmetic Geometric Mean of 1 and 1 / sqrt(2) %% Author: Abhay Jain

-module(agm_calculator). -export([find_agm/0]). -define(TOLERANCE, 0.0000000001).

find_agm() -> 

   A = 1,
   B = 1 / (math:pow(2, 0.5)),
   AGM = agm(A, B),
   io:format("AGM = ~p", [AGM]).

agm (A, B) when abs(A-B) =< ?TOLERANCE -> 

   A;

agm (A, B) -> 

   A1 = (A+B) / 2,
   B1 = math:pow(A*B, 0.5),
   agm(A1, B1).</lang>

Output: <lang Erlang>AGM = 0.8472130848351929</lang>

ERRE

<lang> PROGRAM AGM

! ! for rosettacode.org !

!$DOUBLE

PROCEDURE AGM(A,G->A) 

  LOCAL TA
  REPEAT
     TA=A
     A=(A+G)/2
     G=SQR(TA*G)
  UNTIL A=TA

END PROCEDURE

BEGIN 

  AGM(1.0,1/SQR(2)->A)
  PRINT(A)

END PROGRAM </lang>

F#

Translation of: OCaml

<lang fsharp>let rec agm a g precision = 

   if precision > abs(a - g) then a else
   agm (0.5 * (a + g)) (sqrt (a * g)) precision

printfn "%g" (agm 1. (sqrt(0.5)) 1e-15)</lang> Output 

0.847213

Factor

<lang factor>USING: kernel math math.functions prettyprint ; IN: rosetta-code.arithmetic-geometric-mean

agm ( a g -- a' g' ) 2dup [ + 0.5 * ] 2dip * sqrt ;

1 1 2 sqrt / [ 2dup - 1e-15 > ] [ agm ] while drop .</lang>

Output:
0.8472130847939792

Forth

<lang forth>: agm ( a g -- m ) 

 begin
   fover fover f+ 2e f/
   frot frot f* fsqrt
   fover fover 1e-15 f~
 until
 fdrop ;

1e 2e -0.5e f** agm f. \ 0.847213084793979</lang>

Fortran

A Fortran 77 implementation <lang fortran> function agm(a,b) 

     implicit none
     double precision agm,a,b,eps,c
     parameter(eps=1.0d-15)
  10 c=0.5d0*(a+b)
     b=sqrt(a*b)
     a=c
     if(a-b.gt.eps*a) go to 10
     agm=0.5d0*(a+b)
     end
     program test
     implicit none
     double precision agm
     print*,agm(1.0d0,1.0d0/sqrt(2.0d0))
     end</lang>

FreeBASIC

<lang freebasic>' version 16-09-2015 ' compile with: fbc -s console

Function agm(a As Double, g As Double) As Double 

   Dim As Double t_a 
   
   Do
       t_a = (a + g) / 2
       g = Sqr(a * g)
       Swap a, t_a
   Loop Until a = t_a
   
   Return a

End Function

' ------=< MAIN >=------

Print agm(1, 1 / Sqr(2) )

' empty keyboard buffer While InKey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>

Output:
 0.8472130847939792

Futhark

<lang Futhark> import "futlib/math"

fun agm(a: f64, g: f64): f64 = 

 let eps = 1.0E-16
 loop ((a,g)) = while f64.abs(a-g) > eps do
   ((a+g) / 2.0,
    f64.sqrt (a*g))
 in a

fun main(x: f64, y: f64): f64 = 

 agm(x,y)

</lang>

Go

<lang go>package main

import ( 

   "fmt"
   "math"

)

const ε = 1e-14

func agm(a, g float64) float64 { 

   for math.Abs(a-g) > math.Abs(a)*ε {
       a, g = (a+g)*.5, math.Sqrt(a*g)
   }
   return a

}

func main() { 

   fmt.Println(agm(1, 1/math.Sqrt2))

}</lang>

Output:
0.8472130847939792

Groovy

Translation of: Java

Solution: <lang groovy>double agm (double a, double g) { 

   double an = a, gn = g
   while ((an-gn).abs() >= 10.0**-14) { (an, gn) = [(an+gn)*0.5, (an*gn)**0.5] }
   an

}</lang>

Test: <lang groovy>println "agm(1, 0.5**0.5) = agm(1, ${0.5**0.5}) = ${agm(1, 0.5**0.5)}" assert (0.8472130847939792 - agm(1, 0.5**0.5)).abs() <= 10.0**-14</lang>

Output: 

agm(1, 0.5**0.5) = agm(1, 0.7071067811865476) = 0.8472130847939792

Haskell

<lang haskell>-- Return an approximation to the arithmetic-geometric mean of two numbers. -- The result is considered accurate when two successive approximations are -- sufficiently close, as determined by "eq". agm :: (Floating a) => a -> a -> ((a, a) -> Bool) -> a agm a g eq = snd . head . dropWhile (not . eq) $ iterate step (a, g) 

 where step (a, g) = ((a + g) / 2, sqrt (a * g))

-- Return the relative difference of the pair. We assume that at least one of -- the values is far enough from 0 to not cause problems. relDiff :: (Fractional a) => (a, a) -> a relDiff (x, y) = let n = abs (x - y) 

                    d = ((abs x) + (abs y)) / 2
                in n / d

main = do 

 let equal = (< 0.000000001) . relDiff
 print $ agm 1 (1 / sqrt 2) equal</lang>
Output:
0.8472130847527654

Icon and Unicon

<lang>procedure main(A) 

   a := real(A[1]) | 1.0
   g := real(A[2]) | (1 / 2^0.5)
   epsilon := real(A[3])
   write("agm(",a,",",g,") = ",agm(a,g,epsilon))

end

procedure agm(an, gn, e) 

   /e := 1e-15
   while abs(an-gn) > e do {
      ap := (an+gn)/2.0
      gn := (an*gn)^0.5
      an := ap
      }
   return an

end</lang>

Output: 

->agm
agm(1.0,0.7071067811865475) = 0.8472130847939792
->

J

This one is probably worth not naming, in J, because there are so many interesting variations.

First, the basic approach (with display precision set to 16 digits, which slightly exceeds the accuracy of 64 bit IEEE floating point arithmetic):

<lang j>mean=: +/ % # 

  (mean , */ %:~ #)^:_] 1,%%:2

0.8472130847939792 0.8472130847939791</lang>

This is the limit -- it stops when values are within a small epsilon of previous calculations. We can ask J for unique values (which also means -- unless we specify otherwise -- values within a small epsilon of each other, for floating point values):

<lang j> ~.(mean , */ %:~ #)^:_] 1,%%:2 0.8472130847939792</lang>

Another variation would be to show intermediate values, in the limit process:

<lang j> (mean, */ %:~ #)^:a: 1,%%:2 

                1 0.7071067811865475

0.8535533905932737 0.8408964152537145 0.8472249029234942 0.8472012667468915 0.8472130848351929 0.8472130847527654 0.8472130847939792 0.8472130847939791</lang>

Arbitrary Precision

Another variation would be to use arbitrary precision arithmetic in place of floating point arithmetic.

Borrowing routines from that page, but going with a default of approximately 100 digits of precision:

<lang J>DP=:101

round=: DP&$: : (4 : 0) 

b %~ <.1r2+y*b=. 10x^x

)

sqrt=: DP&$: : (4 : 0) " 0 

assert. 0<:y
%/ <.@%: (2 x: (2*x) round y)*10x^2*x+0>.>.10^.y

)

ln=: DP&$: : (4 : 0) " 0 

assert. 0<y
m=. <.0.5+2^.y
t=. (<:%>:) (x:!.0 y)%2x^m
if. x<-:#":t do. t=. (1+x) round t end.
ln2=. 2*+/1r3 (^%]) 1+2*i.>.0.5*(%3)^.0.5*0.1^x+>.10^.1>.m
lnr=. 2*+/t   (^%]) 1+2*i.>.0.5*(|t)^.0.5*0.1^x
lnr + m * ln2

)

exp=: DP&$: : (4 : 0) " 0 

m=. <.0.5+y%^.2
xm=. x+>.m*10^.2
d=. (x:!.0 y)-m*xm ln 2
if. xm<-:#":d do. d=. xm round d end.
e=. 0.1^xm
n=. e (>i.1:) a (^%!@]) i.>.a^.e [ a=. |y-m*^.2 
(2x^m) * 1++/*/\d%1+i.n

)</lang>

We are also going to want a routine to display numbers with this precision, and we are going to need to manage epsilon manually, and we are going to need an arbitrary root routine:

<lang J>fmt=:[: ;:inv DP&$: : (4 :0)&.> 

 x{.deb (x*2j1)":y

)

root=: ln@] exp@% [

epsilon=: 1r9^DP</lang>

Some example uses:

<lang J> fmt sqrt 2 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572 

  fmt *~sqrt 2

2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 

  fmt epsilon

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000418 

  fmt 2 root 2

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572</lang>

Note that 2 root 2 is considerably slower than sqrt 2. The price of generality. So, while we could define geometric mean generally, a desire for good performance pushes us to use a routine specialized for two numbers:

<lang J>geomean=: */ root~ # geomean2=: [: sqrt */</lang>

A quick test to make sure these can be equivalent:

<lang J> fmt geomean 3 5 3.872983346207416885179265399782399610832921705291590826587573766113483091936979033519287376858673517 

  fmt geomean2 3 5

3.872983346207416885179265399782399610832921705291590826587573766113483091936979033519287376858673517</lang>

Now for our task example:

<lang J> fmt (mean, geomean2)^:(epsilon <&| -/)^:a: 1,%sqrt 2 1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0.707106781186547524400844362104849039284835937688474036588339868995366239231053519425193767163820786 0.853553390593273762200422181052424519642417968844237018294169934497683119615526759712596883581910393 0.840896415253714543031125476233214895040034262356784510813226085974924754953902239814324004199292536 0.847224902923494152615773828642819707341226115600510764553698010236303937284714499763460443890601464 0.847201266746891460403631453693352397963981013612000500823295747923488191871327668107581434542353536 0.847213084835192806509702641168086052652603564606255632688496879079896064578021083935520939216477500 0.847213084752765366704298051779902070392110656059452583317776227659438896688518556753569298762449381 0.847213084793979086607000346473994061522357110332854108003136553369667480633269820344545118989463440 0.847213084793979086605997900490389211440534858586261300461413929971399281619068666682569108141224710 0.847213084793979086606499123482191636481445984459557704232275241670533381126169243513557113565344075 0.847213084793979086606499123482191636481445836194326665888883503648934628542100275932846717790147361 0.847213084793979086606499123482191636481445910326942185060579372659734004834134759723201915677745718 0.847213084793979086606499123482191636481445910326942185060579372659734004834134759723198672311476741 0.847213084793979086606499123482191636481445910326942185060579372659734004834134759723200293994611229 0.847213084793979086606499123482191636481445910326942185060579372659734004834134759723200293994611229</lang>

We could of course extract out only a representative final value, but it's obvious enough, and showing how rapidly this converges is fun.

Java

<lang Java>/* 

* Arithmetic-Geometric Mean of 1 & 1/sqrt(2)
* Brendan Shaklovitz
* 5/29/12
*/

public class ArithmeticGeometricMean { 

   public static double agm(double a, double g) {
       double a1 = a;
       double g1 = g;
       while (Math.abs(a1 - g1) >= 1.0e-14) {
           double arith = (a1 + g1) / 2.0;
           double geom = Math.sqrt(a1 * g1);
           a1 = arith;
           g1 = geom;
       }
       return a1;
   }

   public static void main(String[] args) {
       System.out.println(agm(1.0, 1.0 / Math.sqrt(2.0)));
   }

}</lang>

Output:
0.8472130847939792

JavaScript

ES5

<lang JavaScript>function agm(a0, g0) { 

   var an = (a0 + g0) / 2,
       gn = Math.sqrt(a0 * g0);
   while (Math.abs(an - gn) > tolerance) {
       an = (an + gn) / 2, gn = Math.sqrt(an * gn)
   }
   return an;

}

agm(1, 1 / Math.sqrt(2));</lang>

ES6

<lang JavaScript>(() => { 

   'use strict';

   // ARITHMETIC-GEOMETRIC MEAN

   // agm :: Num a => a -> a -> a
   let agm = (a, g) => {
           let abs = Math.abs,
               sqrt = Math.sqrt;

           return until(
                   m => abs(m.an - m.gn) < tolerance,
                   m => {
                       return {
                           an: (m.an + m.gn) / 2,
                           gn: sqrt(m.an * m.gn)
                       };
                   }, {
                       an: (a + g) / 2,
                       gn: sqrt(a * g)
                   }
               )
               .an;
       },

       // GENERIC

       // until :: (a -> Bool) -> (a -> a) -> a -> a
       until = (p, f, x) => {
           let v = x;
           while (!p(v)) v = f(v);
           return v;
       };



   // TEST

   let tolerance = 0.000001;



   return agm(1, 1 / Math.sqrt(2));

})();</lang>

Output:

<lang JavaScript>0.8472130848351929</lang>

jq

Works with: jq version 1.4

Naive version that assumes tolerance is appropriately specified: <lang jq>def naive_agm(a; g; tolerance): 

 def abs: if . < 0 then -. else . end;
 def _agm:
    # state [an,gn]
    if ((.[0] - .[1])|abs) > tolerance 
    then [add/2, ((.[0] * .[1])|sqrt)] | _agm 
    else .
    end;
 [a, g] | _agm | .[0] ;</lang>

This version avoids an infinite loop if the requested tolerance is too small: <lang jq>def agm(a; g; tolerance): 

 def abs: if . < 0 then -. else . end;
 def _agm:
    # state [an,gn, delta]
    ((.[0] - .[1])|abs) as $delta
    | if $delta == .[2] and $delta < 10e-16 then .
      elif $delta > tolerance
      then [ .[0:2]|add / 2, ((.[0] * .[1])|sqrt), $delta] | _agm
      else . 
      end;
 if tolerance <= 0 then error("specified tolerance must be > 0")
 else [a, g, 0] | _agm | .[0]
 end ;
  1. Example:

agm(1; 1/(2|sqrt); 1e-100)</lang>

Output:
$ jq -n -f Arithmetic-geometric_mean.jq
0.8472130847939792

Julia

<lang Julia>function agm(x::T, y::T, e::Real = 5) where T<:AbstractFloat 

   if x ≤ 0 || y ≤ 0 || e ≤ 0 throw(DomainError("x, y must be strictly positive")) end
   err = e * eps(x)
   g, a = minmax(x, y)
   while err < (a - g)
       a, g = (a + g) / 2, sqrt(a * g)
   end
   return a

end

x = 1.0 y = 1 / √2

println("# Using literal-precision float numbers:") @show agm(x, y) println("# Using half-precision float numbers:") x, y = Float32(x), Float32(y) @show agm(x, y) println("# Using ", precision(BigFloat), "-bit float numbers:") x, y = big(1.0), 1 / √big(2.0) @show agm(x, y)</lang> This version of agm accepts only x and y of matching floating point types. A more permissive version could be created by removing the type parametrization or by creating mixed type versions of the function. The ε for this calculation is given as a positive integer multiple of the machine ε for x.

Output:
# Using literal-precision float numbers:
agm(x, y) = 0.8472130847939792
# Using half-precision float numbers:
agm(x, y) = 0.84721315f0
# Using 256-bit float numbers:
agm(x, y) = 8.472130847939790866064991234821916364814459103269421850605793726597340048341323e-01

Kotlin

<lang scala>// version 1.0.5-2

fun agm(a: Double, g: Double): Double { 

   var aa = a             // mutable 'a'
   var gg = g             // mutable 'g'
   var ta: Double         // temporary variable to hold next iteration of 'aa'
   val epsilon = 1.0e-16  // tolerance for checking if limit has been reached

   while (true) {
       ta = (aa + gg) / 2.0
       if (Math.abs(aa - ta) <= epsilon) return ta
       gg = Math.sqrt(aa * gg)
       aa = ta
   }

}

fun main(args: Array<String>) { 

   println(agm(1.0, 1.0 / Math.sqrt(2.0)))

}</lang>

Output:
0.8472130847939792

LFE

<lang lisp> (defun agm (a g) 

 (agm a g 1.0e-15))

(defun agm (a g tol) 

 (if (=< (- a g) tol)
   a
   (agm (next-a a g)
        (next-g a g)
        tol)))

(defun next-a (a g) 

 (/ (+ a g) 2))

(defun next-g (a g) 

 (math:sqrt (* a g)))

</lang>

Usage: 

> (agm 1 (/ 1 (math:sqrt 2)))
0.8472130847939792

Liberty BASIC

<lang lb> print agm(1, 1/sqr(2)) print using("#.#################",agm(1, 1/sqr(2)))

function agm(a,g) 

   do
       absdiff = abs(a-g)
       an=(a+g)/2
       gn=sqr(a*g)
       a=an
       g=gn
   loop while abs(an-gn)< absdiff
   agm = a

end function

</lang>

LiveCode

<lang LiveCode>function agm aa,g 

   put abs(aa-g) into absdiff 
   put (aa+g)/2 into aan
   put sqrt(aa*g) into gn
   repeat while abs(aan - gn) < absdiff
       put abs(aa-g) into absdiff 
       put (aa+g)/2 into aan
       put sqrt(aa*g) into gn
       put aan into aa
       put gn into g
   end repeat
   return aa

end agm</lang> Example <lang LiveCode>put agm(1, 1/sqrt(2)) -- ouput -- 0.847213</lang>

<lang logo>to about :a :b 

 output and [:a - :b < 1e-15] [:a - :b > -1e-15]

end to agm :arith :geom 

 if about :arith :geom [output :arith]
 output agm (:arith + :geom)/2  sqrt (:arith * :geom)

end

show agm 1 1/sqrt 2 </lang>

Lua

<lang lua>function agm(a, b, tolerance) 

   if not tolerance or tolerance < 1e-15 then
       tolerance = 1e-15
   end
   repeat
       a, b = (a + b) / 2, math.sqrt(a * b)
   until math.abs(a-b) < tolerance
   return a

end

print(string.format("%.15f", agm(1, 1 / math.sqrt(2))))</lang>

Output: 

   0.847213084793979

Maple

Maple provides this function under the name GaussAGM. To compute a floating point approximation, use evalf. <lang Maple> > evalf( GaussAGM( 1, 1 / sqrt( 2 ) ) ); # default precision is 10 digits 

                             0.8472130847

> evalf[100]( GaussAGM( 1, 1 / sqrt( 2 ) ) ); # to 100 digits 0.847213084793979086606499123482191636481445910326942185060579372659\ 

   7340048341347597232002939946112300

</lang> Alternatively, if one or both arguments is already a float, Maple will compute a floating point approximation automatically. <lang Maple> > GaussAGM( 1.0, 1 / sqrt( 2 ) ); 

                             0.8472130847

</lang>

Mathematica

To any arbitrary precision, just increase PrecisionDigits <lang Mathematica>PrecisionDigits = 85; AGMean[a_, b_] := FixedPoint[{ Tr@#/2, Sqrt[Times@@#] }&, N[{a,b}, PrecisionDigits]]〚1〛</lang> 

AGMean[1, 1/Sqrt[2]]
0.8472130847939790866064991234821916364814459103269421850605793726597340048341347597232

MATLAB / Octave

<lang MATLAB>function [a,g]=agm(a,g) %%arithmetic_geometric_mean(a,g) while (1) a0=a; a=(a0+g)/2; g=sqrt(a0*g); if (abs(a0-a) < a*eps) break; end; end; end</lang> 

octave:26> agm(1,1/sqrt(2))
ans =  0.84721

Maxima

<lang maxima>agm(a, b) := %pi/4*(a + b)/elliptic_kc(((a - b)/(a + b))^2)$

agm(1, 1/sqrt(2)), bfloat, fpprec: 85; /* 8.472130847939790866064991234821916364814459103269421850605793726597340048341347597232b-1 */</lang>

МК-61/52

<lang>П1 <-> П0 1 ВП 8 /-/ П2 ИП0 ИП1 - ИП2 - /-/ x<0 31 ИП1 П3 ИП0 ИП1

  • КвКор П1 ИП0 ИП3 + 2 / П0 БП

08 ИП0 С/П</lang>

NetRexx

Translation of: Java

<lang NetRexx>/* NetRexx */ options replace format comments java crossref symbols nobinary

numeric digits 18 parse arg a_ g_ . if a_ = | a_ = '.' then a0 = 1 

                     else a0 = a_

if g_ = | g_ = '.' then g0 = 1 / Math.sqrt(2) 

                     else g0 = g_

say agm(a0, g0)

return

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ method agm(a0, g0) public static returns Rexx 

 a1 = a0
 g1 = g0
 loop while (a1 - g1).abs() >= Math.pow(10, -14)
   temp = (a1 + g1) / 2
   g1 = Math.sqrt(a1 * g1)
   a1 = temp
   end
 return a1 + 0

</lang> Output: 

0.8472130847939792

Nim

<lang nim>import math

proc agm(a, g: float,delta: float = 1.0e-15): float = 

 var 
   aNew: float = 0
   aOld: float = a
   gOld: float = g
 while (abs(aOld - gOld) > delta):
   aNew = 0.5 * (aOld + gOld)
   gOld = sqrt(aOld * gOld)
   aOld = aNew
 result = aOld

echo $agm(1.0,1.0/sqrt(2.0))</lang>

Output:

8.4721308479397917e-01

See first 24 iterations:

<lang nim>from math import sqrt from strutils import parseFloat, formatFloat, ffDecimal

proc agm(x,y: float): tuple[resA,resG: float] = 

 var
   a,g: array[0 .. 23,float]

 a[0] = x
 g[0] = y

 for n in 1 .. 23:
   a[n] = 0.5 * (a[n - 1] + g[n - 1])
   g[n] = sqrt(a[n - 1] * g[n - 1])

 (a[23], g[23])

var t = agm(1, 1/sqrt(2.0))

echo("Result A: " & formatFloat(t.resA, ffDecimal, 24)) echo("Result G: " & formatFloat(t.resG, ffDecimal, 24))</lang>

Oberon-2

Works with: oo2c

<lang oberon2> MODULE Agm; IMPORT 

 Math := LRealMath,
 Out;

CONST 

 epsilon = 1.0E-15;

PROCEDURE Of*(a,g: LONGREAL): LONGREAL; VAR 

 na,ng,og: LONGREAL;

BEGIN 

 na := a; ng := g;
 LOOP
   og := ng;
   ng := Math.sqrt(na * ng); 
   na := (na + og) * 0.5;
   IF na - ng <= epsilon THEN EXIT END
 END;
 RETURN ng;

END Of;

BEGIN 

 Out.LongReal(Of(1,1 / Math.sqrt(2)),0,0);Out.Ln

END Agm. </lang>

Output:
8.4721308479397905E-1

Objeck

Translation of: Java

<lang objeck> class ArithmeticMean { 

 function : Amg(a : Float, g : Float) ~ Nil {
   a1 := a;
   g1 := g;
   while((a1-g1)->Abs() >= Float->Power(10, -14)) {
       tmp := (a1+g1)/2.0;
       g1 := Float->SquareRoot(a1*g1);
       a1 := tmp;
   };
   a1->PrintLine();
 }
 
 function : Main(args : String[]) ~ Nil {
   Amg(1,1/Float->SquareRoot(2));
 }

} </lang>

Output: 

0.847213085

OCaml

<lang ocaml>let rec agm a g tol = 

 if tol > abs_float (a -. g) then a else
 agm (0.5*.(a+.g)) (sqrt (a*.g)) tol

let _ = Printf.printf "%.16f\n" (agm 1.0 (sqrt 0.5) 1e-15)</lang> Output 

0.8472130847939792

Oforth

<lang Oforth>: agm \ a b -- m 

  while( 2dup <> ) [ 2dup + 2 / -rot * sqrt ] drop ;</lang>

Usage : <lang Oforth>1 2 sqrt inv agm</lang>

Output:
0.847213084793979

OOC

<lang ooc> import math // import for sqrt() function

amean: func (x: Double, y: Double) -> Double { 

 (x + y) / 2.

} gmean: func (x: Double, y: Double) -> Double { 

 sqrt(x * y)

} agm: func (a: Double, g: Double) -> Double { 

 while ((a - g) abs() > pow(10, -12)) {
   (a1, g1) := (amean(a, g), gmean(a, g))
   (a, g) = (a1, g1)
 }
 a

}

main: func { 

 "%.16f" printfln(agm(1., sqrt(0.5)))

} </lang> Output 

0.8472130847939792

ooRexx

<lang ooRexx> say agm(1, 1/rxcalcsqrt(2))

routine agm
 use strict arg a, g
 numeric digits 20

 a1 = a
 g1 = g

 loop while abs(a1 - g1) >= 1e-14
     temp = (a1 + g1)/2
     g1 = rxcalcsqrt(a1 * g1)
     a1 = temp
 end
 numeric digits 9
 return a1+0
requires rxmath LIBRARY

</lang>

PARI/GP

Built-in: <lang parigp>agm(1,1/sqrt(2))</lang>

Iteration: <lang parigp>agm2(x,y)=if(x==y,x,agm2((x+y)/2,sqrt(x*y))</lang>

Pascal

Works with: Free_Pascal
Library: GMP

Port of the C example: <lang pascal>Program ArithmeticGeometricMean;

uses 

 gmp;

procedure agm (in1, in2: mpf_t; var out1, out2: mpf_t); begin 

 mpf_add (out1, in1, in2);
 mpf_div_ui (out1, out1, 2);
 mpf_mul (out2, in1, in2);
 mpf_sqrt (out2, out2);

end;

const 

 nl = chr(13)+chr(10);

var 

 x0, y0, resA, resB: mpf_t;
 i: integer;

begin 

 mpf_set_default_prec (65568);

 mpf_init_set_ui (y0, 1);
 mpf_init_set_d (x0, 0.5);
 mpf_sqrt (x0, x0);
 mpf_init (resA);
 mpf_init (resB);

 for i := 0 to 6 do
 begin
   agm(x0, y0, resA, resB);
   agm(resA, resB, x0, y0);
 end;
 mp_printf ('%.20000Ff'+nl, @x0);
 mp_printf ('%.20000Ff'+nl+nl, @y0);

end.</lang> Output is as long as the C example.

Perl

<lang perl>#!/usr/bin/perl -w

my ($a0, $g0, $a1, $g1);

sub agm($$) { 

   $a0 = shift;
   $g0 = shift;
   do { 
       $a1 = ($a0 + $g0)/2; 
       $g1 = sqrt($a0 * $g0); 
       $a0 = ($a1 + $g1)/2; 
       $g0 = sqrt($a1 * $g1); 
   } while ($a0 != $a1); 
   return $a0;

}

print agm(1, 1/sqrt(2))."\n";</lang> Output: 

0.847213084793979

Perl 6

<lang perl6>sub agm( $a is copy, $g is copy ) { 

   ($a, $g) = ($a + $g)/2, sqrt $a * $g until $a ≅ $g;
   return $a;

}

say agm 1, 1/sqrt 2;</lang>

Output:
0.84721308479397917

It's also possible to write it recursively: <lang Perl 6>sub agm( $a, $g ) { 

   $a ≅ $g ?? $a !! agm(|@$_)
       given ($a + $g)/2, sqrt $a * $g;

}

say agm 1, 1/sqrt 2;</lang>

Phix

<lang Phix>function agm(atom a, atom g, atom tolerance=1.0e-15) 

   while abs(a-g)>tolerance do
       {a,g} = {(a + g)/2,sqrt(a*g)}
       printf(1,"%0.15g\n",a)
   end while
   return a

end function ?agm(1,1/sqrt(2)) -- (rounds to 10 d.p.)</lang>

Output:
0.853553390593274
0.847224902923494
0.847213084835193
0.847213084793979
0.8472130848

PHP

<lang php> define('PRECISION', 13);

function agm($a0, $g0, $tolerance = 1e-10) { 

   // the bc extension deals in strings and cannot convert
   // floats in scientific notation by itself - hence
   // this manual conversion to a string
   $limit = number_format($tolerance, PRECISION, '.', );
   $an    = $a0;
   $gn    = $g0;
   do {
       list($an, $gn) = array(
           bcdiv(bcadd($an, $gn), 2), 
           bcsqrt(bcmul($an, $gn)),
       );  
   } while (bccomp(bcsub($an, $gn), $limit) > 0); 

   return $an;

}

bcscale(PRECISION); echo agm(1, 1 / bcsqrt(2)); </lang>

Output:
0.8472130848350

PicoLisp

<lang PicoLisp>(scl 80)

(de agm (A G) 

  (do 7
     (prog1 (/ (+ A G) 2)
        (setq G (sqrt A G)  A @) ) ) )

(round 

  (agm 1.0 (*/ 1.0 1.0 (sqrt 2.0 1.0)))
  70 )</lang>

Output: 

-> "0.8472130847939790866064991234821916364814459103269421850605793726597340"

PL/I

<lang PL/I> arithmetic_geometric_mean: /* 31 August 2012 */ 

  procedure options (main);
  declare (a, g, t) float (18);

  a = 1; g = 1/sqrt(2.0q0);
  put skip list ('The arithmetic-geometric mean of ' || a || ' and ' || g || ':');
  do until (abs(a-g) < 1e-15*a);
     t = (a + g)/2; g = sqrt(a*g);
     a = t;
     put skip data (a, g);
  end;
  put skip list ('The result is:', a);

end arithmetic_geometric_mean; </lang> Results: 

The arithmetic-geometric mean of  1.00000000000000000E+0000 and  7.07106781186547524E-0001: 
A= 8.53553390593273762E-0001                    G= 8.40896415253714543E-0001;
A= 8.47224902923494153E-0001                    G= 8.47201266746891460E-0001;
A= 8.47213084835192807E-0001                    G= 8.47213084752765367E-0001;
A= 8.47213084793979087E-0001                    G= 8.47213084793979087E-0001;
The result is:           8.47213084793979087E-0001

Potion

Input values should be floating point <lang potion>sqrt = (x) : 

  xi = 1
  7 times :
     xi = (xi + x / xi) / 2
  .
  xi

.

agm = (x, y) : 

  7 times :
     a = (x + y) / 2
     g = sqrt(x * y)
     x = a
     y = g
  .
  x

.</lang>

PowerShell

<lang PowerShell> function agm ([Double]$a, [Double]$g) { 

   [Double]$eps = 1E-15
   [Double]$a1 = [Double]$g1 = 0
   while([Math]::Abs($a - $g) -gt $eps) {
       $a1, $g1 = $a, $g
       $a = ($a1 + $g1)/2
       $g = [Math]::Sqrt($a1*$g1)
   }
   [pscustomobject]@{
       a = "$a"
       g = "$g"
   }

} agm 1 (1/[Math]::Sqrt(2)) </lang> Output: 

a                                                      g                                                     
-                                                      -                                                     
0.847213084793979                                      0.847213084793979

Prolog

<lang Prolog> agm(A,G,A) :- abs(A-G) < 1.0e-15, !. agm(A,G,Res) :- A1 is (A+G)/2.0, G1 is sqrt(A*G),!, agm(A1,G1,Res).

?- agm(1,1/sqrt(2),Res). Res = 0.8472130847939792. </lang>

PureBasic

<lang purebasic>Procedure.d AGM(a.d, g.d, ErrLim.d=1e-15) 

 Protected.d ta=a+1, tg
 While ta <> a 
   ta=a: tg=g
   a=(ta+tg)*0.5
   g=Sqr(ta*tg)
 Wend
 ProcedureReturn a

EndProcedure

If OpenConsole() 

 PrintN(StrD(AGM(1, 1/Sqr(2)), 16))
 Input()
 CloseConsole()

EndIf</lang> 

0.8472130847939792

Python

The calculation generates two new values from two existing values which is the classic example for the use of assignment to a list of values in the one statement, so ensuring an gn are only calculated from an-1 gn-1.

Basic Version

<lang python>from math import sqrt

def agm(a0, g0, tolerance=1e-10): 

   """
   Calculating the arithmetic-geometric mean of two numbers a0, g0.

   tolerance     the tolerance for the converged 
                 value of the arithmetic-geometric mean
                 (default value = 1e-10)
   """
   an, gn = (a0 + g0) / 2.0, sqrt(a0 * g0)
   while abs(an - gn) > tolerance:
       an, gn = (an + gn) / 2.0, sqrt(an * gn)
   return an

print agm(1, 1 / sqrt(2))</lang>

Output:
 0.847213084835

Multi-Precision Version

<lang python>from decimal import Decimal, getcontext

def agm(a, g, tolerance=Decimal("1e-65")): 

   while True:
       a, g = (a + g) / 2, (a * g).sqrt()
       if abs(a - g) < tolerance:
           return a

getcontext().prec = 70 print agm(Decimal(1), 1 / Decimal(2).sqrt())</lang>

Output:
0.847213084793979086606499123482191636481445910326942185060579372659734

All the digits shown are correct.

R

<lang r>arithmeticMean <- function(a, b) { (a + b)/2 } geometricMean <- function(a, b) { sqrt(a * b) }

arithmeticGeometricMean <- function(a, b) { 

 rel_error <- abs(a - b) / pmax(a, b) 
 if (all(rel_error < .Machine$double.eps, na.rm=TRUE)) {
   agm <- a
   return(data.frame(agm, rel_error));
 }
 Recall(arithmeticMean(a, b), geometricMean(a, b))

}

agm <- arithmeticGeometricMean(1, 1/sqrt(2)) print(format(agm, digits=16))</lang>

Output:
                 agm             rel_error
1 0.8472130847939792 1.310441309927519e-16

This function also works on vectors a and b (following the spirit of R): <lang r>a <- c(1, 1, 1) b <- c(1/sqrt(2), 1/sqrt(3), 1/2) agm <- arithmeticGeometricMean(a, b) print(format(agm, digits=16))</lang>

Output:
                 agm             rel_error
1 0.8472130847939792 1.310441309927519e-16
2 0.7741882646460426 0.000000000000000e+00
3 0.7283955155234534 0.000000000000000e+00

Racket

This version uses Racket's normal numbers: <lang racket>

  1. lang racket

(define (agm a g [ε 1e-15]) 

 (if (<= (- a g) ε)
     a
     (agm (/ (+ a g) 2) (sqrt (* a g)) ε)))

(agm 1 (/ 1 (sqrt 2))) </lang> Output: 

0.8472130847939792

This alternative version uses arbitrary precision floats: <lang racket>

  1. lang racket

(require math/bigfloat) (bf-precision 200) (bfagm 1.bf (bf/ (bfsqrt 2.bf))) </lang> Output: 

(bf #e0.84721308479397908660649912348219163648144591032694218506057918)

Raven

<lang Raven>define agm use $a, $g, $errlim 

   # $errlim $g $a "%d %g %d\n" print
   $a 1.0  +   as $t
   repeat $a 1.0 * $g - abs -15 exp10 $a *  >   while
       $a $g + 2 /   as $t
       $a $g * sqrt  as $g
       $t as $a
       $g $a $t  "t: %g a: %g g: %g\n" print
   $a


16 1 2 sqrt / 1 agm "agm: %.15g\n" print</lang>

Output:
t: 0.853553 a: 0.853553 g: 0.840896
t: 0.847225 a: 0.847225 g: 0.847201
t: 0.847213 a: 0.847213 g: 0.847213
t: 0.847213 a: 0.847213 g: 0.847213
agm: 0.847213084793979

REXX

Also, this version of the AGM REXX program has three   short circuits   within it for an equality case and for two zero cases.

REXX supports arbitrary precision, so the default digits can be changed if desired. <lang rexx>/*REXX program calculates the AGM (arithmetic─geometric mean) of two (real) numbers. */ parse arg a b digs . /*obtain optional numbers from the C.L.*/ if digs== | digs=="," then digs=110 /*No DIGS specified? Then use default.*/ numeric digits digs /*REXX will use lots of decimal digits.*/ if a== | a=="," then a=1 /*No A specified? Then use the default*/ if b== | b=="," then b=1 / sqrt(2) /* " B " " " " " */ call AGM a,b /*invoke the AGM function. */ say '1st # =' a /*display the A value. */ say '2nd # =' b /* " " B " */ say ' AGM =' agm(a, b) /* " " AGM " */ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ agm: procedure: parse arg x,y; if x=y then return x /*is this an equality case?*/ 

                                 if y=0  then return 0      /*is   Y  equal to zero ?  */
                                 if x=0  then return y / 2  /* "   X    "    "   "     */
     d=digits();  numeric digits d+5            /*add 5 more digs to ensure convergence*/
     tiny='1e-' || (digits() - 1)               /*construct a pretty tiny REXX number. */
     ox=x + 1
                       do #=1  while ox\=x & abs(ox)>tiny;  ox=x;          oy=y
                                                             x=(ox+oy)/2;   y=sqrt(ox*oy)
                       end   /*#*/
     numeric digits d;       return x / 1       /*restore digs, normalize X to new digs*/

/*──────────────────────────────────────────────────────────────────────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); m.=9; numeric form; h=d+6 

     numeric digits; parse value format(x,2,1,,0) 'E0'  with  g 'E' _ .;  g=g *.5'e'_ % 2
       do j=0  while h>9;      m.j=h;               h=h % 2  + 1;  end /*j*/
       do k=j+5  to 0  by -1;  numeric digits m.k;  g=(g+x/g)*.5;  end /*k*/;    return g</lang>
output   when using the default input:
1st # = 1
2nd # = 0.70710678118654752440084436210484903928483593768847403658833986899536623923105351942519376716382078636750692312
  AGM = 0.84721308479397908660649912348219163648144591032694218506057937265973400483413475972320029399461122994212228563

Ring

<lang ring> decimals(9) see agm(1, 1/sqrt(2)) + nl see agm(1,1/pow(2,0.5)) + nl

func agm agm,g 

      while agm
            an  = (agm + g)/2
            gn  = sqrt(agm*g)
            if fabs(agm-g) <= fabs(an-gn) exit ok
            agm = an
            g   = gn
      end
      return gn

</lang>

Ruby

Flt Version

The thing to note about this implementation is that it uses the Flt library for high-precision math. This lets you adapt context (including precision and epsilon) to a ridiculous-in-real-life degree. <lang ruby># The flt package (http://flt.rubyforge.org/) is useful for high-precision floating-point math.

  1. It lets us control 'context' of numbers, individually or collectively -- including precision
  2. (which adjusts the context's value of epsilon accordingly).

require 'flt' include Flt

BinNum.Context.precision = 512 # default 53 (bits)

def agm(a,g) 

 new_a = BinNum a
 new_g = BinNum g
 while new_a - new_g > new_a.class.Context.epsilon do
   old_g = new_g
   new_g = (new_a * new_g).sqrt
   new_a = (old_g + new_a) * 0.5
 end
 new_g

end

puts agm(1, 1 / BinNum(2).sqrt)</lang>

Output:
0.84721308479397908660649912348219163648144591032694218506057937265973400483413475972320029399461122994212228562523341096309796266583087105969971363598338426

Adjusting the precision setting (at about line 9) will of course affect this. :-)

BigDecimal Version

Ruby has a BigDecimal class in standard library <lang ruby>require 'bigdecimal'

PRECISION = 100 EPSILON = 0.1 ** (PRECISION/2) BigDecimal::limit(PRECISION)

def agm(a,g) 

 while a - g > EPSILON
   a, g = (a+g)/2, (a*g).sqrt(PRECISION)
 end
 [a, g]

end

a = BigDecimal(1) g = 1 / BigDecimal(2).sqrt(PRECISION) puts agm(a, g)</lang>

Output:
0.847213084793979086606499123482191636481445910326942185060579372659734004834134759723201915677745718E0
0.8472130847939790866064991234821916364814459103269421850605793726597340048341347597231986723114767413E0

Run BASIC

<lang runbasic>print agm(1, 1/sqr(2)) print agm(1,1/2^.5) print using("#.############################",agm(1, 1/sqr(2)))

function agm(agm,g) 

while agm
  an  = (agm + g)/2
  gn  = sqr(agm*g)
  if abs(agm-g) <= abs(an-gn) then exit while
  agm = an
  g   = gn
wend

end function</lang>Output: 

0.847213085
0.847213085
0.8472130847939791165772005376

Rust

<lang rust>// Accepts two command line arguments // cargo run --name agm arg1 arg2

fn main () { 

   let mut args = std::env::args();

   let x = args.nth(1).expect("First argument not specified.").parse::<f32>().unwrap();
   let y = args.next().expect("Second argument not specified.").parse::<f32>().unwrap();

   let result = agm(x,y);
   println!("The arithmetic-geometric mean is {}", result);

}

fn agm (x: f32, y: f32) -> f32 { 

   let e: f32 = 0.000001;
   let mut a = x;
   let mut g = y;
   let mut a1: f32;
   let mut g1: f32;

   if a * g < 0f32 { panic!("The arithmetric-geometric mean is undefined for numbers less than zero!"); }
   else {
       loop {
           a1 = (a + g) / 2.;
           g1 = (a * g).sqrt();
           a = a1;
           g = g1;
           if (a - g).abs() < e {  return a; }
       }
   }

}</lang>

Output:

Output of running with arguments 1, 0.70710678: 

The arithmetic-geometric mean is 1.456791

Scala

<lang scala> 

 def agm(a: Double, g: Double, eps: Double): Double = {
   if (math.abs(a - g) < eps) (a + g) / 2
   else agm((a + g) / 2, math.sqrt(a * g), eps)
 }

 agm(1, math.sqrt(2)/2, 1e-15)

</lang>

Scheme

<lang scheme> (define agm 

 (case-lambda
   ((a0 g0) ; call again with default value for tolerance
    (agm a0 g0 1e-8))
   ((a0 g0 tolerance) ; called with three arguments
    (do ((a a0 (* (+ a g) 1/2))
         (g g0 (sqrt (* a g))))
      ((< (abs (- a g)) tolerance) a)))))

(display (agm 1 (/ 1 (sqrt 2)))) (newline) </lang>

Output:
0.8472130848351929

Seed7

<lang seed7>$ include "seed7_05.s7i"; 

 include "float.s7i";
 include "math.s7i";

const func float: agm (in var float: a, in var float: g) is func 

 result
   var float: agm is 0.0;
 local
   const float: iota is 1.0E-7;
   var float: a1 is 0.0;
   var float: g1 is 0.0;
 begin
   if a * g < 0.0 then
     raise RANGE_ERROR;
   else
     while abs(a - g) > iota do
       a1 := (a + g) / 2.0;
       g1 := sqrt(a * g);
       a := a1;
       g := g1;
     end while;
     agm := a;
   end if;
 end func;

const proc: main is func 

 begin
   writeln(agm(1.0, 2.0) digits 6);
   writeln(agm(1.0, 1.0 / sqrt(2.0)) digits 6);
 end func;</lang>
Output:
1.456791
0.847213

SequenceL

<lang sequencel>import <Utilities/Math.sl>;

agm(a, g) := 

   let
       iota := 1.0e-15;
       arithmeticMean := 0.5 * (a + g);
       geometricMean := sqrt(a * g);
   in
       a when abs(a-g) < iota
   else
       agm(arithmeticMean, geometricMean);

main := agm(1.0, 1.0 / sqrt(2));</lang>

Output:
0.847213

Sidef

<lang ruby>func agm(a, g) { 

   loop {
       var (a1, g1) = ((a+g)/2, sqrt(a*g))
       [a1,g1] == [a,g] && return a
       (a, g) = (a1, g1)
   }

}

say agm(1, 1/sqrt(2))</lang>

Output:
0.8472130847939790866064991234821916364814

Sinclair ZX81 BASIC

Translation of: COBOL

Works with 1k of RAM.

The specification calls for a function. Sadly that is not available to us, so this program uses a subroutine: pass the arguments in the global variables A and G, and the result will be returned in AGM. The performance is quite acceptable. Note that the subroutine clobbers A and G, so you should save them if you want to use them again.

Better precision than this is not easily obtainable on the ZX81, unfortunately. <lang basic> 10 LET A=1 

20 LET G=1/SQR 2
30 GOSUB 100
40 PRINT AGM
50 STOP

100 LET A0=A 110 LET A=(A+G)/2 120 LET G=SQR (A0*G) 130 IF ABS(A-G)>.00000001 THEN GOTO 100 140 LET AGM=A 150 RETURN</lang>

Output:
0.84721309

Stata

<lang stata>mata

real scalar agm(real scalar a, real scalar b) { real scalar c do { c=0.5*(a+b) b=sqrt(a*b) a=c } while (a-b>1e-15*a) return(0.5*(a+b)) }

agm(1,1/sqrt(2)) end</lang>

Output:
.8472130848

Swift

<lang Swift>import Darwin

enum AGRError : Error { case undefined }

func agm(_ a: Double, _ g: Double, _ iota: Double = 1e-8) throws -> Double { var a = a var g = g var a1: Double = 0 var g1: Double = 0

guard a * g >= 0 else { throw AGRError.undefined }

while abs(a - g) > iota { a1 = (a + g) / 2 g1 = sqrt(a * g) a = a1 g = g1 }

return a }

do { try print(agm(1, 1 / sqrt(2))) } catch { print("agr is undefined when a * g < 0") }</lang>

Output:
0.847213084835193

Tcl

The tricky thing about this implementation is that despite the finite precision available to IEEE doubles (which Tcl uses in its implementation of floating point arithmetic, in common with many other languages) the sequence of values does not quite converge to a single value; it gets to within a ULP and then errors prevent it from getting closer. This means that an additional termination condition is required: once a value does not change (hence the old_b variable) we have got as close as we can. Note also that we are using exact equality with floating point; this is reasonable because this is a rapidly converging sequence (it only takes 4 iterations in this case). <lang tcl>proc agm {a b} { 

   set old_b [expr {$b<0?inf:-inf}]
   while {$a != $b && $b != $old_b} {

set old_b $b lassign [list [expr {0.5*($a+$b)}] [expr {sqrt($a*$b)}]] a b 

   }
   return $a

}

puts [agm 1 [expr 1/sqrt(2)]]</lang> Output: 

0.8472130847939792

TI-83 BASIC

<lang ti83b>1→A:1/sqrt(2)→G While abs(A-G)>e-15 (A+G)/2→B sqrt(AG)→G:B→A End A</lang>

Output:
.8472130848

UNIX Shell

Works with: ksh93

ksh is one of the few unix shells that can do floating point arithmetic (bash does not). <lang bash>function agm { 

   float a=$1 g=$2 eps=${3:-1e-11} tmp
   while (( abs(a-g) > eps )); do
       print "debug: a=$a\tg=$g"
       tmp=$(( (a+g)/2.0 ))
       g=$(( sqrt(a*g) ))
       a=$tmp
   done
   echo $a

}

agm $((1/sqrt(2))) 1</lang>

Output:
debug: a=0.7071067812	g=1
debug: a=0.8535533906	g=0.8408964153
debug: a=0.8472249029	g=0.8472012668
debug: a=0.8472130848	g=0.8472130847
debug: a=0.8472130848	g=0.8472130848
debug: a=0.8472130848	g=0.8472130848
0.8472130848

You can get a more approximate convergence by changing the while condition to compare the numbers as strings: change <lang bash>while (( abs(a-g) > eps ))</lang> to <lang bash>while $a != $g </lang>

VBScript

Translation of: BBC BASIC

<lang vb> Function agm(a,g) Do Until a = tmp_a tmp_a = a a = (a + g)/2 g = Sqr(tmp_a * g) Loop agm = a End Function

WScript.Echo agm(1,1/Sqr(2)) </lang>

Output:
0.847213084793979

XPL0

<lang XPL0>include c:\cxpl\codesi; real A, A1, G; [Format(0, 16); A:= 1.0; G:= 1.0/sqrt(2.0); repeat A1:= (A+G)/2.0; G:= sqrt(A*G); A:= A1; RlOut(0, A); RlOut(0, G); RlOut(0, A-G); CrLf(0); until A=G; ]</lang>

Output: 

 8.5355339059327400E-001 8.4089641525371500E-001 1.2656975339559100E-002
 8.4722490292349400E-001 8.4720126674689100E-001 2.3636176602726000E-005
 8.4721308483519300E-001 8.4721308475276500E-001 8.2427509262572600E-011
 8.4721308479397900E-001 8.4721308479397900E-001 0.0000000000000000E+000

zkl

Translation of: XPL0

<lang zkl>a:=1.0; g:=1.0/(2.0).sqrt(); while(not a.closeTo(g,1.0e-15)){ 

  a1:=(a+g)/2.0; g=(a*g).sqrt(); a=a1; 
  println(a,"  ",g," ",a-g);

}</lang>

Output:
0.853553  0.840896 0.012657
0.847225  0.847201 2.36362e-05
0.847213  0.847213 8.24275e-11
0.847213  0.847213 1.11022e-16

Or, using tail recursion <lang zkl>fcn(a=1.0, g=1.0/(2.0).sqrt()){ println(a," ",g," ",a-g); 

  if(a.closeTo(g,1.0e-15)) return(a) else return(self.fcn((a+g)/2.0, (a*g).sqrt()));

}()</lang>

Output:
1 0.707107 0.292893
0.853553 0.840896 0.012657
0.847225 0.847201 2.36362e-05
0.847213 0.847213 8.24275e-11
0.847213 0.847213 1.11022e-16

ZX Spectrum Basic

Translation of: ERRE

<lang zxbasic>10 LET a=1: LET g=1/SQR 2 20 LET ta=a 30 LET a=(a+g)/2 40 LET g=SQR (ta*g) 50 IF a<ta THEN GO TO 20 60 PRINT a </lang>

Output:
0.84721309
Retrieved from "https://rosettacode.org/wiki/Arithmetic-geometric_mean?oldid=212246"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information