

 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)



Task

Write a function to compute the arithmetic-geometric mean of two numbers.



The arithmetic-geometric mean of two numbers can be (usefully) denoted as 




a
g
m

(
a
,
g
)


{\displaystyle \mathrm {agm} (a,g)}

, and is equal to the limit of the sequence:

 




a

0


=
a
;


g

0


=
g


{\displaystyle a_{0}=a;\qquad g_{0}=g}


 




a

n
+
1


=



1
2



(

a

n


+

g

n


)
;


g

n
+
1


=



a

n



g

n




.


{\displaystyle a_{n+1}={\tfrac {1}{2}}(a_{n}+g_{n});\quad g_{n+1}={\sqrt {a_{n}g_{n}}}.}


Since the limit of 




a

n


−

g

n




{\displaystyle a_{n}-g_{n}}

 tends (rapidly) to zero with iterations, this is an efficient method.

Demonstrate the function by calculating:






a
g
m

(
1
,
1

/



2


)


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




Also see

   mathworld.wolfram.com/Arithmetic-Geometric Mean


MODULE AGM;
FROM EXCEPTIONS IMPORT AllocateSource,ExceptionSource,GetMessage,RAISE;
FROM LongConv IMPORT ValueReal;
FROM LongMath IMPORT sqrt;
FROM LongStr IMPORT RealToStr;
FROM Terminal IMPORT ReadChar,Write,WriteString,WriteLn;
 
VAR
    TextWinExSrc : ExceptionSource;
 
PROCEDURE ReadReal() : LONGREAL;
VAR
    buffer : ARRAY[0..63] OF CHAR;
    i : CARDINAL;
    c : CHAR;
BEGIN
    i := 0;
 
    LOOP
        c := ReadChar();
        IF ((c >= '0') AND (c <= '9')) OR (c = '.') THEN
            buffer[i] := c;
            Write(c);
            INC(i)
        ELSE
            WriteLn;
            EXIT
        END
    END;
 
    buffer[i] := 0C;
    RETURN ValueReal(buffer)
END ReadReal;
 
PROCEDURE WriteReal(r : LONGREAL);
VAR
    buffer : ARRAY[0..63] OF CHAR;
BEGIN
    RealToStr(r, buffer);
    WriteString(buffer)
END WriteReal;
 
PROCEDURE AGM(a,g : LONGREAL) : LONGREAL;
CONST iota = 1.0E-16;
VAR a1, g1 : LONGREAL;
BEGIN
    IF a * g < 0.0 THEN
        RAISE(TextWinExSrc, 0, "arithmetic-geometric mean undefined when x*y<0")
    END;
 
    WHILE ABS(a - g) > iota DO
        a1 := (a + g) / 2.0;
        g1 := sqrt(a * g);
 
        a := a1;
        g := g1
    END;
 
    RETURN a
END AGM;
 
VAR
    x, y, z: LONGREAL;
BEGIN
    WriteString("Enter two numbers: ");
    x := ReadReal();
    y := ReadReal();
    WriteReal(AGM(x, y));
    WriteLn
END AGM.