Arithmetic-geometric mean

Arithmetic-geometric mean is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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)

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 $\mathrm {agm} (a,g)$, and is equal to the limit of the sequence:

$a_{0}=a;\qquad g_{0}=g$

$a_{n+1}={\tfrac {1}{2}}(a_{n}+g_{n});\quad g_{n+1}={\sqrt {a_{n}g_{n}}}.$

Since the limit of $a_{n}-g_{n}$ tends (rapidly) to zero with iterations, this is an efficient method.

Demonstrate the function by calculating:

$\mathrm {agm} (1,1/{\sqrt {2}})$

C++

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

 Nigel_Galloway
 February 7th., 2011.

• /

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.

J

This one is 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>

Another variation would be to use arbitrary precision arithmetic in place of floating point arithmetic. (out of time, maybe later)

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