Arithmetic-geometric mean

From Rosetta Code
Jump to: navigation, search
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 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 angn tends (rapidly) to zero with iterations, this is an efficient method.

Demonstrate the function by calculating:

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

Contents

[edit] 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;
Output:
0.84721308479397909

[edit] AutoHotkey

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

Output:

0.847213084793979


[edit] 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)
}
 

Output

0.8472130847939792

[edit] BBC BASIC

      *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
 

Produces this output:

0.8472130847939792

[edit] 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)
Output:
.84721308479397908659

[edit] C

#include<math.h>
#include<stdio.h>
#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;
}
 


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

[edit] C#

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));
}
}
}
}

Output:

0.847213084835193

[edit] C

/*Arithmetic Geometric Mean of 1 and 1/sqrt(2)
 
Nigel_Galloway
February 7th., 2012.
*/

 
#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;
}

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.

[edit] Common 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)))
 
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

[edit] 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)));
}
Output:
0.8472130847939790866

All the digits shown are exact.

[edit] 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).

Output:

AGM = 0.8472130848351929


[edit] F#

Translation of: OCaml
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)

Output

0.847213

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

[edit] Fortran

A Fortran 77 implementation

      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

[edit] 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))
}
Output:
0.8472130847939792

[edit] 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
Output:
0.8472130847527654

[edit] Icon and Unicon

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

Output:

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

[edit] 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):

mean=: +/ % #
(mean , */ %:~ #)^:_] 1,%%:2
0.8472130847939792 0.8472130847939791

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):

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

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

   (mean, */ %:~ #)^:a: 1,%%:2
1 0.7071067811865475
0.8535533905932737 0.8408964152537145
0.8472249029234942 0.8472012667468915
0.8472130848351929 0.8472130847527654
0.8472130847939792 0.8472130847939791

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

[edit] 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));
}
}
Output:
(0.8472130847939792)

[edit] 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));

[edit] Liberty BASIC

 
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
 
 

[edit]

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
 

[edit] Maple

Maple provides this function under the name GaussAGM. To compute a floating point approximation, use evalf.

 
> 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
 

Alternatively, if one or both arguments is already a float, Maple will compute a floating point approximation automatically.

 
> GaussAGM( 1.0, 1 / sqrt( 2 ) );
0.8472130847
 

[edit] Mathematica

To any arbitrary precision, just increase PrecisionDigits

PrecisionDigits = 85;
AGMean[a_, b_] := FixedPoint[{ Tr@#/2, Sqrt[Times@@#] }&, N[{a,b}, PrecisionDigits]]〚1〛
AGMean[1, 1/Sqrt[2]]
0.8472130847939790866064991234821916364814459103269421850605793726597340048341347597232

[edit] MATLAB / Octave

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
octave:26> agm(1,1/sqrt(2))
ans =  0.84721

[edit] 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 */

[edit] МК-61/52

П1	<->	П0	1	ВП	8	/-/	П2	ИП0	ИП1
- ИП2 - /-/ x<0 31 ИП1 П3 ИП0 ИП1
* КвКор П1 ИП0 ИП3 + 2 / П0 БП
08 ИП0 С/П

[edit] NetRexx

Translation of: Java
/* 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
 

Output:

0.8472130847939792

[edit] Nimrod

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

Output:

8.4721308479397917e-01

[edit] Objeck

Translation of: Java
 
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));
}
}
 

Output:

0.847213085

[edit] 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)

Output

0.8472130847939792

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

[edit] PARI/GP

Built-in:

agm(1,1/sqrt(2))

Iteration:

agm2(x,y)=if(x==y,x,agm2((x+y)/2,sqrt(x*y))

[edit] Pascal

Works with: Free_Pascal
Library: GMP

Port of the C example:

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.

Output is as long as the C example.

[edit] 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";

Output:

0.847213084793979

[edit] Perl 6

 
sub agm( $a is copy, $g is copy ) {
loop {
given ($a + $g)/2, sqrt $a * $g {
return $a if @$_ ~~ ($a, $g);
($a, $g) = @$_;
}
}
}
 
say agm 1, 1/sqrt 2;
 
 
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.

It's also possible to write it recursively:

 
sub agm( $a, $g ) {
@$_ ~~ ($a, $g) ?? $a !! agm(|@$_)
given ($a + $g)/2, sqrt $a * $g;
}
 
say agm 1, 1/sqrt 2;

[edit] 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));
 
Output:
0.8472130848350

[edit] 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 )

Output:

-> "0.8472130847939790866064991234821916364814459103269421850605793726597340"

[edit] 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;
 

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 

[edit] 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
0.8472130847939792

[edit] Python

[edit] Basic Version

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))
Output:
 0.847213084835

[edit] Multi-Precision Version

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())
Output:
0.847213084793979086606499123482191636481445910326942185060579372659734

All the digits shown are correct.

[edit] 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))
Output:
                 agm             rel_error
1 0.8472130847939792 1.310441309927519e-16

This function also works on vectors a and b (following the spirit of 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))
Output:
                 agm             rel_error
1 0.8472130847939792 1.310441309927519e-16
2 0.7741882646460426 0.000000000000000e+00
3 0.7283955155234534 0.000000000000000e+00

[edit] Racket

This version uses Racket's normal numbers:

 
#lang racket
(define (agm a g [ε 1e-15])
(if (<= (- a g) ε)
a
(agm (/ (+ a g) 2) (sqrt (* a g)) ε)))
 
(agm 1 (/ 1 (sqrt 2)))
 

Output:

0.8472130847939792

This alternative version uses arbitrary precision floats:

 
#lang racket
(require math/bigfloat)
(bf-precision 200)
(bfagm 1.bf (bf/ (bfsqrt 2.bf)))
 

Output:

(bf #e0.84721308479397908660649912348219163648144591032694218506057918)

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

[edit] 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 for two zero cases.
REXX supports arbitrary precision, so the digits can be increased if desired.

/*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 /*stick a fork in it, we're done.*/
/*────────────────────────────AGM subroutine────────────────────────────*/
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 /*while*/
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

output when using the default input:

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

[edit] 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.

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

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


[edit] Run BASIC

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
Output:
0.847213085
0.847213085
0.8472130847939791165772005376

[edit] 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)
 

[edit] 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;
Output:
1.456791
0.847213

[edit] 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).

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)]]

Output:

0.8472130847939792

[edit] 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;
]

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

[edit] zkl

Translation of: XPL0
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);
}
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

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()));
}()
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
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox