Arithmetic-geometric mean: Difference between revisions

=={{header|11l}}==

{{trans|Python}}

V an = (a0 + g0) / 2.0

V gn = sqrt(a0 * g0)

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

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

USING AGM,R13

SAVEAREA B STM-SAVEAREA(R15)

 

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

 

with: n

{{libheader|Action! Tool Kit}}

{{libheader|Action! Real Math}}

 

PROC Agm(REAL POINTER a0,g0,result)

 

=={{header|Ada}}==

 

procedure Arith_Geom_Mean is

 

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

BEGIN

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

 

=={{header|APL}}==

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

1 agd ÷2*÷2

By functional composition:

 

 

property tolerance : 1.0E-5

 

=={{header|AutoHotkey}}==

While abs(a-g) > tolerance

{

 

=={{header|AWK}}==

BEGIN {

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

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

{{works with|QBasic}}

END

 

 

==={{header|BASIC256}}===

end

 

 

==={{header|Commodore BASIC}}===

20 G = 1/SQR(2)

30 GOSUB 100

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

{{works with|BBC BASIC for Windows}}

@% = &1010

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

 

==={{header|GW-BASIC}}===

20 G = 1!/SQR(2!)

30 FOR I=1 TO 20 'twenty iterations is plenty

 

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

110 DEF AGM(A,G)

120 DO

==={{header|True BASIC}}===

{{works with|QBasic}}

DO

LET ta = (a + g) / 2

 

=={{header|bc}}==

* numbers x and y.

* Result will have d digits after the decimal point.

 

=={{header|BQN}}==

(|𝕨-𝕩) ≤ 1e¯15? 𝕨;

(0.5×𝕨+𝕩) 𝕊 √𝕨×𝕩

=={{header|C}}==

===Basic===

#include<stdio.h>

#include<stdlib.h>

 

===GMP===

 

Nigel_Galloway

 

=={{header|C sharp|C#}}==

{

using System;

Note that the last 5 digits are spurious, as ''maximumRelativeDifference'' was only specified to be 1e-5. Using 1e-11 instead will give the result 0.847213084793979, which is as far as ''double'' can take it.

===Using Decimal Type===

 

class Program {

{{Libheader|System.Numerics}}

Even though the System.Numerics library directly supports only '''BigInteger''' (and not big rationals or big floating point numbers), it can be coerced into making this calculation. One just has to keep track of the decimal place and multiply by a very large constant.

using static System.Console;

using BI = System.Numerics.BigInteger;

 

=={{header|C++}}==

#include<bits/stdc++.h>

using namespace std;

 

=={{header|Clojure}}==

(:gen-class))

 

 

=={{header|COBOL}}==

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

DATA DIVISION.

 

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

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

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

 

=={{header|D}}==

 

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

{{libheader| System.SysUtils}}

{{Trans|C#}}

program geometric_mean;

 

=={{header|EchoLisp}}==

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

(lib 'math)

 

=={{header|Elixir}}==

 

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

def mean(a,g,tol) do

 

=={{header|Erlang}}==

%% Author: Abhay Jain

 

agm(A1, B1).</syntaxhighlight>

Output:

 

=={{header|ERRE}}==

PROGRAM AGM

 

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

{{trans|OCaml}}

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

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

 

=={{header|Factor}}==

IN: rosetta-code.arithmetic-geometric-mean

 

 

=={{header|Forth}}==

begin

fover fover f+ 2e f/

=={{header|Fortran}}==

A '''Fortran 77''' implementation

implicit none

double precision agm,a,b,eps,c

 

=={{header|FreeBASIC}}==

' compile with: fbc -s console

 

{{incorrect|Futhark|Futhark's syntax has changed, so this example will not compile}}

 

import "futlib/math"

 

 

=={{header|Go}}==

 

import (

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

 

Test:

assert (0.8472130847939792 - agm(1, 0.5**0.5)).abs() <= 10.0**-14</syntaxhighlight>

 

 

=={{header|Haskell}}==

-- The result is considered accurate when two successive approximations are

-- sufficiently close, as determined by "eq".

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

 

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

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

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

0.8472130847939792 0.8472130847939791</syntaxhighlight>

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

 

0.8472130847939792</syntaxhighlight>

 

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

 

1 0.7071067811865475

0.8535533905932737 0.8408964152537145

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

 

 

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

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

)

Some example uses:

 

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572

fmt *~sqrt 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:

 

geomean2=: [: sqrt */</syntaxhighlight>

 

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

=={{header|Java}}==

 

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

* Brendan Shaklovitz

 

===ES5===

var an = (a0 + g0) / 2,

gn = Math.sqrt(a0 * g0);

 

===ES6===

'use strict';

 

 

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

[a, g] | _agm | .[0] ;</syntaxhighlight>

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

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

def _agm:

=={{header|Julia}}==

{{works with|Julia|1.2}}

(x ≤ 0 || y ≤ 0 || e ≤ 0) && throw(DomainError("x, y must be strictly positive"))

g, a = minmax(x, y)

=={{header|Klingphix}}==

{{trans|Oforth}}

 

:agm [ over over + 2 / rot rot * sqrt ] [ over over tostr swap tostr # ] while drop ;

" " input</syntaxhighlight>

{{trans|F#}}

 

:agm %a %g %p !p !g !a

 

=={{header|Kotlin}}==

 

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

=={{header|LFE}}==

 

(defun agm (a g)

(agm a g 1.0e-15))

 

=={{header|Liberty BASIC}}==

print agm(1, 1/sqr(2))

print using("#.#################",agm(1, 1/sqr(2)))

 

=={{header|LiveCode}}==

put abs(aa-g) into absdiff

put (aa+g)/2 into aan

end agm</syntaxhighlight>

Example

-- ouput

-- 0.847213</syntaxhighlight>

 

=={{header|LLVM}}==

; LLVM does not provide a way to print values, so the alternative would be

; to just load the string into memory, and that would be boring.

 

=={{header|Logo}}==

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

end

=={{header|Lua}}==

 

if not tolerance or tolerance < 1e-15 then

tolerance = 1e-15

 

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

Module Checkit {

Function Agm {

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

</syntaxhighlight>

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

=={{header|Mathematica}}/{{header|Wolfram Language}}==

To any arbitrary precision, just increase PrecisionDigits

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

 

 

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

%%arithmetic_geometric_mean(a,g)

while (1)

 

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

{{trans|C}}

FROM EXCEPTIONS IMPORT AllocateSource,ExceptionSource,GetMessage,RAISE;

FROM LongConv IMPORT ValueReal;

=={{header|NetRexx}}==

{{trans|Java}}

options replace format comments java crossref symbols nobinary

 

 

=={{header|NewLISP}}==

(define (a-next a g) (mul 0.5 (add a g)))

 

 

=={{header|Nim}}==

 

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

See first 24 iterations:

 

from strutils import parseFloat, formatFloat, ffDecimal

 

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

{{works with|oo2c}}

MODULE Agm;

IMPORT

=={{header|Objeck}}==

{{trans|Java}}

class ArithmeticMean {

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

 

=={{header|OCaml}}==

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

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

=={{header|Oforth}}==

 

while( 2dup <> ) [ 2dup + 2 / -rot * sqrt ] drop ;</syntaxhighlight>

 

Usage :

 

{{out}}

 

=={{header|OOC}}==

import math // import for sqrt() function

 

 

=={{header|ooRexx}}==

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

 

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

Built-in:

 

Iteration:

 

=={{header|Pascal}}==

{{libheader|GMP}}

Port of the C example:

 

uses

 

=={{header|Perl}}==

 

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

 

=={{header|Phix}}==

<span style="color: #008080;">function</span> <span style="color: #000000;">agm</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">g</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">tolerance</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1.0e-15</span><span style="color: #0000FF;">)</span>

<span style="color: #008080;">while</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">-</span><span style="color: #000000;">g</span><span style="color: #0000FF;">)></span><span style="color: #000000;">tolerance</span> <span style="color: #008080;">do</span>

 

=={{header|Phixmonti}}==

 

1.0e-15 var tolerance

 

=={{header|PHP}}==

define('PRECISION', 13);

 

 

=={{header|Picat}}==

println(agm(1.0, 1/sqrt(2))).

 

 

=={{header|PicoLisp}}==

 

(de agm (A G)

 

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

arithmetic_geometric_mean: /* 31 August 2012 */

procedure options (main);

=={{header|Potion}}==

Input values should be floating point

xi = 1

7 times :

 

=={{header|PowerShell}}==

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

[Double]$eps = 1E-15

 

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

 

=={{header|PureBasic}}==

Protected.d ta=a+1, tg

While ta <> a

 

===Basic Version===

 

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

<pre> 0.847213084835</pre>

===Multi-Precision Version===

 

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

=={{header|Quackery}}==

 

 

[ temp put

 

=={{header|R}}==

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

 

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)

=={{header|Racket}}==

This version uses Racket's normal numbers:

#lang racket

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

 

This alternative version uses arbitrary precision floats:

#lang racket

(require math/bigfloat)

=={{header|Raku}}==

(formerly Perl 6)

($a, $g) = ($a + $g)/2, sqrt $a * $g until $a ≅ $g;

return $a;

 

It's also possible to write it recursively:

$a ≅ $g ?? $a !! agm(|@$_)

given ($a + $g)/2, sqrt $a * $g;

 

=={{header|Raven}}==

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

$a 1.0 + as $t

 

=={{header|Relation}}==

function agm(x,y)

set a = x

 

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

if digs=='' | digs=="," then digs= 120 /*No DIGS specified? Then use default.*/

 

=={{header|Ring}}==

decimals(9)

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

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

===BigDecimal Version===

Ruby has a BigDecimal class in standard library

 

PRECISION = 100

 

=={{header|Run BASIC}}==

print agm(1,1/2^.5)

print using("#.############################",agm(1, 1/sqr(2)))

=={{header|Rust}}==

 

// cargo run --name agm arg1 arg2

 

 

=={{header|Scala}}==

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

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

=={{header|Scheme}}==

 

(define agm

(case-lambda

 

=={{header|Seed7}}==

include "float.s7i";

include "math.s7i";

 

=={{header|SequenceL}}==

 

agm(a, g) :=

 

=={{header|Sidef}}==

loop {

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

 

Better precision than this is not easily obtainable on the ZX81, unfortunately.

20 LET G=1/SQR 2

30 GOSUB 100

{{works with|Smalltalk/X}}

That is simply a copy/paste of the already existing agm method in the Number class:

"return the arithmetic-geometric mean agm(x, y)

of the receiver (x) and the argument, y.

^ ai</syntaxhighlight>

 

Transcript showCR: ( (1/2) agm:(1/6) ).

Transcript showCR: (1 agm:(1 / 2 sqrt)).</syntaxhighlight>

{{works with|oracle|11.2 and higher}}

The solution uses recursive WITH clause (aka recursive CTE, recursive query, recursive factored subquery). Some, perhaps many, but not all SQL dialects support recursive WITH clause. The solution below was written and tested in Oracle SQL - Oracle has supported recursive WITH clause since version 11.2.

rec (rn, a, g, diff) as (

select 1, 1, 1/sqrt(2), 1 - 1/sqrt(2)

 

=={{header|Standard ML}}==

fun agm(a, g) = let

fun agm'(a, g, eps) =

 

=={{header|Stata}}==

 

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

 

=={{header|Swift}}==

 

enum AGRError : Error {

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

 

=={{header|TI-83 BASIC}}==

While abs(A-G)>e-15

(A+G)/2→B

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

0.8472130848</pre>

 

 

=={{header|VBA}}==

Do While Abs(a - g) > tolerance

tmp = a

=={{header|VBScript}}==

{{trans|BBC BASIC}}

Function agm(a,g)

Do Until a = tmp_a

{{trans|C#}}

===Double, Decimal Versions===

Imports System.Console

 

{{trans|C#}}

{{Libheader|System.Numerics}}

Imports System.Console

Imports BI = System.Numerics.BigInteger

 

=={{header|Vlang}}==

const ep = 1e-14

}</syntaxhighlight>

Using standard math module

import math

 

=={{header|Wren}}==

{{trans|Go}}

 

var agm = Fn.new { |a, g|

 

=={{header|XPL0}}==

real A, A1, G;

[Format(0, 16);

=={{header|zkl}}==

{{trans|XPL0}}

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

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

</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()));

}()</syntaxhighlight>

=={{header|ZX Spectrum Basic}}==

{{trans|ERRE}}

20 LET ta=a

30 LET a=(a+g)/2