Arithmetic-geometric mean: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎{{header|ooRexx}}: Add example for ooRexx)
m (→‎{{header|REXX}}: changed the name of an internal subroutine for a subroutine to be consistant with other REXX examples. -- ~~~~)
Line 405: Line 405:
/*───────────────────────────────SQRT subroutine────────────────────────*/
/*───────────────────────────────SQRT subroutine────────────────────────*/
sqrt: procedure; parse arg x;if x=0 then return 0;d=digits();numeric digits 11
sqrt: procedure; parse arg x;if x=0 then return 0;d=digits();numeric digits 11
g=$sqguess(); do j=0 while p>9; m.j=p; p=p%2+1; end
g=.sqrtGuess(); do j=0 while p>9; m.j=p; p=p%2+1; end
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k
g=.5*(g+x/g); end; numeric digits d; return g/1
g=.5*(g+x/g); end; numeric digits d; return g/1


$sqguess: numeric form scientific;m.=11; p=d+d%4+2
.sqrtGuess: numeric form scientific;m.=11; p=d+d%4+2
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2</lang>
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2</lang>
'''output''' when using the default input:
'''output''' when using the default input:

Revision as of 18:24, 11 June 2012

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)

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:

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

C++

<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.

D

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

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

   do {
       //(a, g) = tuple((a + g) / 2.0, sqrt(a * g));
       immutable ag = tuple((a + g) / 2.0, sqrt(a * g));
       a = ag[0];
       g = ag[1];
   } while (feqrel(a, g) < bitPrecision);
   return a;

}

void main() {

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

}</lang>

Output:
0.8472130847939790866

All the digits shown are exact.

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

The referenced Mathworld page mentions that AGM is meaningful on the complex plane as well.

It certainly has. It has been called The Mind of God (perhaps beyond the scope of this task!). There is a little more to it than adding +0i to the real solution. In Real Maths sqrt 4 is 2(This task assumes this). In fantasy maths sqrt 4 is 2 and -2. Then the next iteration takes the sqrt of a negative number, and maths goes complex. The result in Complex maths is the set of all these arithmetric geometric means.--Nigel Galloway 14:20, 23 April 2012 (UTC)

<lang go>package main

import (

   "fmt"
   "math"
   "math/cmplx"

)

const ε = 1e-14

func agm(a, g complex128) complex128 {

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

}

func main() {

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

}</lang>

Output:
(0.8472130847939792+0i)

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)

Java

<lang Java>/*

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

public class ArithmeticMean {

       public static void agm (double a, double g){
               double a1 = a;
               double g1 = g;
               while (Math.abs(a1-g1) >= Math.pow(10, -14)){
                       double aTemp = (a1+g1)/2.0;
                       g1 = Math.sqrt(a1*g1);
                       a1 = aTemp;
               }
               System.out.println(a1);
       }
       public static void main(String[] args){
               agm(1,1/Math.sqrt(2));
       }

}</lang>

Output:
(0.8472130847939792)

Mathematica

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

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

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>

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>

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 6

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

   sub iter ($old) {
       my $new := [ 0.5 * [+](@$old), sqrt [*](@$old) ];
       last if $new ~~ $old;
       $new;
   }
   ([$a,$g], &iter ... 0)[*-1][0];

}

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

Output:
0.84721308479397917

Obviously the "fixed point" detector here is relying on the floating-point representation running out of bits, or this algorithm would not terminate before using up all memory.

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"

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

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.

REXX

The REXX language doesn't have a SQRT function, so one was included here (RYO).
Also, this version of AGM has three short circuits within it for an equality case and and two zero cases.
REXX supports arbitrary precision, so that can be increased if desired. <lang rexx>/*REXX program calculates AGM (arithmetric-geometric mean) of 2 numbers.*/

parse arg a b digs . /*obtain numbers from the command line.*/ if digs== then digs=100 /*no DIGS specified? Then use default.*/ numeric digits digs /*Now, REXX will use lots of digits. */ if a== then a=1 /*no A specified? Then use default. */ if b== then b=1/sqrt(2) /*no B specified? " " " */ say '1st # =' a say '2nd # =' b say ' AGM =' agm(a,b)/1 /*divide by 1; goes from 105-->100 digs*/ say ' AGM =' agm(a,b)/1 /*dividing by 1 normalizes the REXX num*/ exit /*───────────────────────────────AGM (arithmetric─geometric mean) sub.──*/ agm: procedure: parse arg x,y; if x=y then return x /*equality case.*/ if y=0 then return 0; if x=0 then return .5*y /*two "0" cases.*/ numeric digits digits()+5 /*add 5 more digs to ensure convergence*/ !='1e-' || (digits()-1); _x=x+1

               do while _x\=x & abs(_x)>!;   _x=x;   _y=y;   x=(_x+_y)*.5
               y=sqrt(_x*_y)
               end

return x /*───────────────────────────────SQRT subroutine────────────────────────*/ sqrt: procedure; parse arg x;if x=0 then return 0;d=digits();numeric digits 11

                g=.sqrtGuess();  do j=0 while p>9;  m.j=p;  p=p%2+1;  end
                do k=j+5 to 0 by -1;  if m.k>11 then numeric digits m.k
                g=.5*(g+x/g);  end;  numeric digits d;  return g/1

.sqrtGuess: numeric form scientific;m.=11; p=d+d%4+2

         parse value format(x,2,1,,0) 'E0' with g 'E' _ .;  return g*.5'E'_%2</lang>

output when using the default input:

1st # = 1
2nd # = 0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207862
  AGM = 0.8472130847939790866064991234821916364814459103269421850605793726597340048341347597232002939946112299

Ruby

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>

  1. The flt package (http://flt.rubyforge.org/) is useful for high-precision floating-point math.
  2. It lets us control 'context' of numbers, individually or collectively -- including precision
  3. (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. :-)


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

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>

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

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