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>

* &nbsp; [http://mathworld.wolfram.com/Arithmetic-GeometricMean.html mathworld.wolfram.com/Arithmetic-Geometric Mean]

<br><br>

 

=={{header|360 Assembly}}==

For maximum compatibility, this program uses only the basic instruction set.

USING AGM,R13

SAVEAREA B STM-SAVEAREA(R15)

LTORG

YREGS

{{out}}

<pre>

 

=={{header|8th}}==

 

with: n

;with

bye

{{out}}

<pre>

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

</pre>

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

 

Printing out the difference between the means at each iteration nicely demonstrates the quadratic convergence.

BEGIN

PROC agm = (LONG REAL x, y) LONG REAL :

ELIF x + y = LONG 0.0 THEN LONG 0.0 CO Edge cases CO

ELSE

FI

END;

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

END

Output:<pre>+1.707106781186547524400844362e +0 +2.928932188134524755991556379e -1

+2.928932188134524755991556379e -1 +1.265697533955921916929670477e -2

 

=={{header|APL}}==

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

1 agd ÷2*÷2

Output: <pre>0.8472130848</pre>

 

By functional composition:

 

 

property tolerance : 1.0E-5

end script

end if

{{Out}}

<pre>0.847213084835</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>

 

=={{header|BASIC}}==

 

==={{header|BBC BASIC}}===

{{works with|BBC BASIC for Windows}}

@% = &1010

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

UNTIL a = ta

= a

 

==={{header|IS-BASIC}}===

110 DEF AGM(A,G)

120 DO

150 LOOP UNTIL A=TA

160 LET AGM=A

 

=={{header|bc}}==

* numbers x and y.

* Result will have d digits after the decimal point.

 

scale = 20

 

{{Out}}

<pre>.84721308479397908659</pre>

 

=={{header|C}}==

===Basic===

#include<stdio.h>

#include<stdlib.h>

return 0;

}

 

Original output:

 

===GMP===

 

Nigel_Galloway

 

return 0;

 

The first couple of iterations produces:

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.

 

{

using System;

}

}

Output:

<pre>0.847213084835193</pre>

 

=={{header|Clojure}}==

(:gen-class))

 

 

(println (agm one isqrt2))

{{Output}}

<pre>

 

=={{header|COBOL}}==

PROGRAM-ID. ARITHMETIC-GEOMETRIC-MEAN-PROG.

DATA DIVISION.

COMPUTE G = FUNCTION SQRT(G).

SUBTRACT A FROM G GIVING DIFF.

{{out}}

<pre>0.8472130847939792</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}}==

 

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

void main() @safe {

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

{{out}}

<pre>0.8472130847939790866</pre>

All the digits shown are exact.

 

=={{header|EchoLisp}}==

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

(lib 'math)

 

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

→ 0.8472130847939792

 

=={{header|Elixir}}==

 

def mean(a,g,tol) when abs(a-g) <= tol, do: a

def mean(a,g,tol) do

end

 

 

{{out}}

 

=={{header|Erlang}}==

%% Author: Abhay Jain

 

A1 = (A+B) / 2,

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

Output:

 

=={{header|ERRE}}==

PROGRAM AGM

 

PRINT(A)

END PROGRAM

 

=={{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|Factor}}==

IN: rosetta-code.arithmetic-geometric-mean

 

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

 

{{out}}

<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|Futhark}}==

 

import "futlib/math"

 

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

agm(x,y)

 

=={{header|Go}}==

 

import (

func main() {

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

{{out}}

<pre>

{{trans|Java}}

Solution:

double an = a, gn = g

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

an

 

Test:

 

Output:

 

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

 

=== Arbitrary Precision ===

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

 

 

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

n=. e (>i.1:) a (^%!@]) i.>.a^.e [ a=. |y-m*^.2

(2x^m) * 1++/*/\d%1+i.n

 

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

 

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

)

root=: ln@] exp@% [

 

 

Some example uses:

 

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572

fmt *~sqrt 2

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000418

fmt 2 root 2

 

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

 

 

A quick test to make sure these can be equivalent:

 

3.872983346207416885179265399782399610832921705291590826587573766113483091936979033519287376858673517

fmt geomean2 3 5

 

Now for our task example:

 

1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0.707106781186547524400844362104849039284835937688474036588339868995366239231053519425193767163820786

0.853553390593273762200422181052424519642417968844237018294169934497683119615526759712596883581910393 0.840896415253714543031125476233214895040034262356784510813226085974924754953902239814324004199292536

0.847213084793979086606499123482191636481445984459557704232275241670533381126169243513557113565344075 0.847213084793979086606499123482191636481445836194326665888883503648934628542100275932846717790147361

0.847213084793979086606499123482191636481445910326942185060579372659734004834134759723201915677745718 0.847213084793979086606499123482191636481445910326942185060579372659734004834134759723198672311476741

 

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

=={{header|Java}}==

 

* Arithmetic-Geometric Mean of 1 & 1/sqrt(2)

* Brendan Shaklovitz

System.out.println(agm(1.0, 1.0 / Math.sqrt(2.0)));

}

{{out}}

<pre>0.8472130847939792</pre>

 

===ES5===

var an = (a0 + g0) / 2,

gn = Math.sqrt(a0 * g0);

}

 

 

===ES6===

'use strict';

 

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

 

 

{{Out}}

 

=={{header|jq}}==

{{works with|jq|1.4}}

Naive version that assumes tolerance is appropriately specified:

def abs: if . < 0 then -. else . end;

def _agm:

else .

end;

This version avoids an infinite loop if the requested tolerance is too small:

def abs: if . < 0 then -. else . end;

def _agm:

# Example:

{{Out}}

$ jq -n -f Arithmetic-geometric_mean.jq

 

=={{header|Julia}}==

g, a = minmax(x, y)

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

end

end

println("# Using literal-precision float numbers:")

@show agm(x, y)

println("# Using half-precision float numbers:")

x, y = Float32(x), Float32(y)

@show agm(x, y)

println("# Using ", precision(BigFloat), "-bit float numbers:")

x, y = big(1.0), 1 / √big(2.0)

 

{{out}}

# Using 256-bit float numbers:

agm(x, y) = 8.472130847939790866064991234821916364814459103269421850605793726597340048341323e-01</pre>

 

=={{header|Kotlin}}==

 

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

fun main(args: Array<String>) {

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

 

{{out}}

0.8472130847939792

</pre>

 

=={{header|LFE}}==

 

(defun agm (a g)

(agm a g 1.0e-15))

(defun next-g (a g)

(math:sqrt (* a g)))

 

Usage:

0.8472130847939792

</pre>

 

=={{header|LiveCode}}==

put abs(aa-g) into absdiff

put (aa+g)/2 into aan

end repeat

return aa

Example

-- ouput

 

=={{header|Logo}}==

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

end

 

show agm 1 1/sqrt 2

 

=={{header|Lua}}==

 

if not tolerance or tolerance < 1e-15 then

tolerance = 1e-15

end

 

 

'''Output:'''

 

=={{header|M2000 Interpreter}}==

Module Checkit {

Function Agm {

}

Checkit

 

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

0.8472130847939792

</pre>

 

=={{header|Nim}}==

 

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

 

 

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

 

=={{header|Oberon-2}}==

{{works with|oo2c}}

MODULE Agm;

IMPORT

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

END Agm.

{{Out}}

<pre>

=={{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|Oforth}}==

 

 

Usage :

 

{{out}}

 

=={{header|OOC}}==

import math // import for sqrt() function

 

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

}

Output

<pre>0.8472130847939792</pre>

 

=={{header|ooRexx}}==

say agm(1, 1/rxcalcsqrt(2,16))

 

return a1+0

 

{{out}}

<pre>0.8472130847939791968</pre>

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

 

=={{header|Phix}}==

{{out}}

<pre>

0.8472130848

</pre>

 

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

=={{header|Potion}}==

Input values should be floating point

xi = 1

7 times :

.

x

 

=={{header|PowerShell}}==

function agm ([Double]$a, [Double]$g) {

[Double]$eps = 1E-15

}

agm 1 (1/[Math]::Sqrt(2))

<b>Output:</b>

<pre>

 

=={{header|Prolog}}==

agm(A,G,A) :- abs(A-G) < 1.0e-15, !.

agm(A,G,Res) :- A1 is (A+G)/2.0, G1 is sqrt(A*G),!, agm(A1,G1,Res).

?- agm(1,1/sqrt(2),Res).

Res = 0.8472130847939792.

 

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

 

REXX supports arbitrary precision, so the default digits can be changed if desired.

parse arg a b digs . /*obtain optional numbers from the C.L.*/

numeric digits digs /*REXX will use lots of decimal digits.*/

call AGM a,b /*invoke the AGM function. */

say '1st # =' a /*display the A value. */

say '2nd # =' b /* " " B " */

say ' AGM =' agm(a, b) /* " " AGM " */

/*──────────────────────────────────────────────────────────────────────────────────────*/

d= digits() /*obtain the current decimal digits. */

numeric digits d + 5 /*add 5 more digs to ensure convergence*/

ox= x + 1 /*ensure that the old X ¬= new X. */

do while ox\=x & abs(ox)>tiny /*compute until the old X ≡ new X. */

x= (ox + oy) * .5 /*compute " new " " X. */

y= sqrt(ox * oy) /* " " " " " Y. */

numeric digits; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g *.5'e'_ % 2

do j=0 while h>9; m.j=h; h=h % 2 + 1; end /*j*/

{{out|output|text=&nbsp; when using the default input:}}

<pre>

1st # = 1

</pre>

 

=={{header|Ring}}==

decimals(9)

see agm(1, 1/sqrt(2)) + nl

end

return gn

 

=={{header|Ruby}}==

===Flt Version===

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

end

 

{{out}}

<pre>0.84721308479397908660649912348219163648144591032694218506057937265973400483413475972320029399461122994212228562523341096309796266583087105969971363598338426</pre>

===BigDecimal Version===

Ruby has a BigDecimal class in standard library

 

PRECISION = 100

a = BigDecimal(1)

g = 1 / BigDecimal(2).sqrt(PRECISION)

{{out}}

<pre>

0.8472130847939790866064991234821916364814459103269421850605793726597340048341347597231986723114767413E0

</pre>

 

=={{header|Rust}}==

// cargo run --name agm arg1 arg2

 

}

}

 

{{out}}

 

=={{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|Scheme}}==

 

(define agm

(case-lambda

 

(display (agm 1 (/ 1 (sqrt 2)))) (newline)

 

{{out}}

 

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

 

=={{header|SequenceL}}==

 

agm(a, g) :=

agm(arithmeticMean, geometricMean);

 

 

{{out}}

 

=={{header|Sidef}}==

loop {

var (a1, g1) = ((a+g)/2, sqrt(a*g))

}

 

{{out}}

<pre>0.8472130847939790866064991234821916364814</pre>

 

 

 

{{out}}

 

=={{header|Stata}}==

 

real scalar agm(real scalar a, real scalar b) {

 

agm(1,1/sqrt(2))

{{out}}

<pre>.8472130848</pre>

 

=={{header|Swift}}==

 

enum AGRError : Error {

} catch {

print("agr is undefined when a * g < 0")

{{out}}

<pre>0.847213084835193</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>

 

{{out}}

 

=={{header|UNIX Shell}}==

{{works with|ksh93}}

ksh is one of the few unix shells that can do floating point arithmetic (bash does not).

float a=$1 g=$2 eps=${3:-1e-11} tmp

while (( abs(a-g) > eps )); do

}

 

 

{{output}}

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

0.847213 0.847213 1.11022e-16

</pre>