Arithmetic-geometric mean: Difference between revisions

Write a function to compute the [[wp:Arithmetic-geometric mean|arithmetic-geometric mean]] of two numbers.

The arithmetic-geometric mean of two numbers can be (usefully) denoted as <math>\mathrm{agm}(a,g)</math>, and is equal to the limit of the sequence:

: <math>a_0 = a; \qquad g_0 = g</math>

Demonstrate the function by calculating:

:<math>\mathrm{agm}(1,1/\sqrt{2})</math>

 

=={{header|Ada}}==

 

procedure Arith_Geom_Mean is

begin

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

 

Output:<pre>0.84721308479397909</pre>

 

=={{header|AutoHotkey}}==

While abs(a-g) > tolerance

{

}

SetFormat, FloatFast, 0.15

Output:

<pre>0.847213084793979</pre>

 

=={{header|AWK}}==

BEGIN {

printf "%.16g\n", agm(1.0,sqrt(0.5))

return (x<0 ? -x : x)

}

Output

<pre>0.8472130847939792</pre>

 

{{works with|BBC BASIC for Windows}}

@% = &1010

PRINT FNagm(1, 1/SQR(2))

UNTIL a = ta

= a

 

=={{header|bc}}==

* numbers x and y.

* Result will have d digits after the decimal point.

 

scale = 20

 

{{Out}}

<pre>.84721308479397908659</pre>

 

=={{header|C}}==

#include<stdio.h>

#include<stdlib.h>

return 0;

}

 

<pre>

Enter two numbers: 1.0 2.0

The arithmetic-geometric mean is 1.456791

</pre>

 

{

using System;

}

}

Output:

<pre>0.847213084835193</pre>

 

 

 

 

}

 

 

 

 

 

 

<pre>

</pre>

 

 

=={{header|Common Lisp}}==

(loop for a = a0 then (* (+ a g) 5d-1)

and g = g0 then (sqrt (* a g))

until (< (abs (- a g)) tolerance)

finally (return a)))

 

{{out}}

CL-USER> (agm 1d0 (/ 1d0 (sqrt 2d0)) 1d-12)

0.8472130847939792d0</pre>

=={{header|D}}==

 

do {

} while (feqrel(a, g) < bitPrecision);

return a;

}

 

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

{{out}}

<pre>0.8472130847939790866</pre>

All the digits shown are exact.

 

=={{header|Erlang}}==

%% Author: Abhay Jain

 

A1 = (A+B) / 2,

B1 = math:pow(A*B, 0.5),

Output:

 

 

=={{header|F_Sharp|F#}}==

{{trans|OCaml}}

if precision > abs(a - g) then a else

agm (0.5 * (a + g)) (sqrt (a * g)) precision

 

Output

<pre>0.847213</pre>

 

=={{header|Forth}}==

begin

fover fover f+ 2e f/

fdrop ;

 

 

=={{header|Fortran}}==

A '''Fortran 77''' implementation

implicit none

double precision agm,a,b,eps,c

double precision agm

print*,agm(1.0d0,1.0d0/sqrt(2.0d0))

 

=={{header|Go}}==

 

import (

func main() {

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

{{out}}

<pre>

0.8472130847939792

</pre>

 

=={{header|Haskell}}==

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

 

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

 

main = do

let equal = (< 0.000000001) . relDiff

{{out}}

<pre>0.8472130847527654</pre>

=={{header|Icon}} and {{header|Unicon}}==

 

a := real(A[1]) | 1.0

g := real(A[2]) | (1 / 2^0.5)

}

return an

 

Output:

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

 

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

 

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

 

 

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

 

1 0.7071067811865475

0.8535533905932737 0.8408964152537145

0.8472249029234942 0.8472012667468915

0.8472130848351929 0.8472130847527654

 

 

 

 

 

 

 

 

}

{{out}}

 

=={{header|JavaScript}}==

}

}

 

 

 

 

=={{header|Logo}}==

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

end

 

show agm 1 1/sqrt 2

 

=={{header|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

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

 

To any arbitrary precision, just increase PrecisionDigits

 

<pre>AGMean[1, 1/Sqrt[2]]

 

=={{header|MATLAB}} / {{header|Octave}}==

%%arithmetic_geometric_mean(a,g)

while (1)

if (abs(a0-a) < a*eps) break; end;

end;

<pre>octave:26> agm(1,1/sqrt(2))

ans = 0.84721

 

=={{header|Maxima}}==

 

agm(1, 1/sqrt(2)), bfloat, fpprec: 85;

 

=={{header|МК-61/52}}==

- ИП2 - /-/ x<0 31 ИП1 П3 ИП0 ИП1

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

 

=={{header|NetRexx}}==

{{trans|Java}}

options replace format comments java crossref symbols nobinary

 

end

return a1 + 0

'''Output:'''

<pre>

</pre>

 

 

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

result = aOld

 

 

Output:<br/>

See first 24 iterations:

 

from strutils import parseFloat, formatFloat, ffDecimal

 

(a[23], g[23])

 

 

echo("Result A: " & formatFloat(t.resA, ffDecimal, 24))

 

=={{header|Objeck}}==

{{trans|Java}}

class ArithmeticMean {

function : Amg(a : Float, g : Float) ~ Nil {

}

}

 

Output:

 

=={{header|OCaml}}==

if tol > abs_float (a -. g) then a else

agm (0.5*.(a+.g)) (sqrt (a*.g)) tol

 

Output

<pre>0.8472130847939792</pre>

 

=={{header|OOC}}==

import math // import for sqrt() function

 

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

}

Output

<pre>0.8472130847939792</pre>

 

=={{header|ooRexx}}==

 

::routine agm

loop while abs(a1 - g1) >= 1e-14

temp = (a1 + g1)/2

a1 = temp

end

return a1+0

 

 

=={{header|PARI/GP}}==

Built-in:

 

Iteration:

 

=={{header|Pascal}}==

{{libheader|GMP}}

Port of the C example:

 

uses

mp_printf ('%.20000Ff'+nl, @x0);

mp_printf ('%.20000Ff'+nl+nl, @y0);

Output is as long as the C example.

 

=={{header|Perl}}==

 

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

}

 

Output:

<pre>0.847213084793979</pre>

 

{{out}}

 

 

 

=={{header|PHP}}==

define('PRECISION', 13);

 

bcscale(PRECISION);

echo agm(1, 1 / bcsqrt(2));

{{out}}

<pre>

0.8472130848350

</pre>

 

=={{header|PicoLisp}}==

 

(de agm (A G)

(round

(agm 1.0 (*/ 1.0 1.0 (sqrt 2.0 1.0)))

Output:

<pre>-> "0.8472130847939790866064991234821916364814459103269421850605793726597340"</pre>

 

=={{header|PL/I}}==

arithmetic_geometric_mean: /* 31 August 2012 */

procedure options (main);

put skip list ('The result is:', a);

end arithmetic_geometric_mean;

Results:

<pre>

</pre>

 

 

 

 

=={{header|Python}}==

===Basic Version===

 

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

return an

 

{{out}}

<pre> 0.847213084835</pre>

===Multi-Precision Version===

 

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

 

getcontext().prec = 70

{{out}}

<pre>0.847213084793979086606499123482191636481445910326942185060579372659734</pre>

All the digits shown are correct.

 

=={{header|R}}==

geometricMean <- function(a, b) { sqrt(a * b) }

 

 

agm <- arithmeticGeometricMean(1, 1/sqrt(2))

{{out}}

<pre> agm rel_error

1 0.8472130847939792 1.310441309927519e-16</pre>

This function also works on vectors a and b (following the spirit of R):

b <- c(1/sqrt(2), 1/sqrt(3), 1/2)

agm <- arithmeticGeometricMean(a, b)

{{out}}

<pre> agm rel_error

=={{header|Racket}}==

This version uses Racket's normal numbers:

#lang racket

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

 

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

Output:

<pre>

 

This alternative version uses arbitrary precision floats:

#lang racket

(require math/bigfloat)

(bf-precision 200)

(bfagm 1.bf (bf/ (bfsqrt 2.bf)))

Output:

<pre>

(bf #e0.84721308479397908660649912348219163648144591032694218506057918)

</pre>

 

=={{header|Raven}}==

# $errlim $g $a "%d %g %d\n" print

$a 1.0 + as $t

 

 

{{out}}

<pre>t: 0.853553 a: 0.853553 g: 0.840896

t: 0.847213 a: 0.847213 g: 0.847213

agm: 0.847213084793979</pre>

 

=={{header|REXX}}==

 

end /*while*/

 

1st # = 1

</pre>

 

=={{header|Ruby}}==

The thing to note about this implementation is that it uses the [http://flt.rubyforge.org/ Flt] library for high-precision math. This lets you adapt context (including precision and epsilon) to a ridiculous-in-real-life degree.

# It lets us control 'context' of numbers, individually or collectively -- including precision

# (which adjusts the context's value of epsilon accordingly).

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

 

new_a = BinNum a

new_g = BinNum g

end

 

{{out}}

<pre>0.84721308479397908660649912348219163648144591032694218506057937265973400483413475972320029399461122994212228562523341096309796266583087105969971363598338426</pre>

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

 

 

 

 

=={{header|Scala}}==

def agm(a: Double, g: Double, eps: Double): Double = {

if (math.abs(a - g) < eps) (a + g) / 2

 

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

 

=={{header|Seed7}}==

include "float.s7i";

include "math.s7i";

writeln(agm(1.0, 2.0) digits 6);

writeln(agm(1.0, 1.0 / sqrt(2.0)) digits 6);

 

{{out}}

<pre>

1.456791

0.847213

</pre>

 

=={{header|Sidef}}==

}

}

 

{{out}}

<pre>0.8472130847939790866064991234821916364814</pre>

 

=={{header|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 <code>old_b</code> 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).

set old_b [expr {$b<0?inf:-inf}]

while {$a != $b && $b != $old_b} {

}

 

Output:

<pre>0.8472130847939792</pre>

 

=={{header|XPL0}}==

real A, A1, G;

[Format(0, 16);

RlOut(0, A); RlOut(0, G); RlOut(0, A-G); CrLf(0);

until A=G;

 

Output:

=={{header|zkl}}==

{{trans|XPL0}}

while(not a.closeTo(g,1.0e-15)){

a1:=(a+g)/2.0; g=(a*g).sqrt(); a=a1;

println(a," ",g," ",a-g);

{{out}}

<pre>

</pre>

Or, using tail recursion

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

{{out}}

<pre>