

 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


AGM      CSECT  
         USING  AGM,R13
SAVEAREA B      STM-SAVEAREA(R15)
         DC     17F'0'
         DC     CL8'AGM'
STM      STM    R14,R12,12(R13)
         ST     R13,4(R15)
         ST     R15,8(R13)
         LR     R13,R15
         ZAP    A,K                a=1
         ZAP    PWL8,K
         MP     PWL8,K
         DP     PWL8,=P'2'
         ZAP    PWL8,PWL8(7)
         BAL    R14,SQRT
         ZAP    G,PWL8             g=sqrt(1/2)
WHILE1   EQU    *                  while a!=g
         ZAP    PWL8,A
         SP     PWL8,G
         CP     PWL8,=P'0'         (a-g)!=0
         BE     EWHILE1
         ZAP    PWL8,A
         AP     PWL8,G
         DP     PWL8,=P'2'
         ZAP    AN,PWL8(7)         an=(a+g)/2
         ZAP    PWL8,A
         MP     PWL8,G
         BAL    R14,SQRT
         ZAP    G,PWL8             g=sqrt(a*g)
         ZAP    A,AN               a=an
         B      WHILE1
EWHILE1  EQU    *
         ZAP    PWL8,A
         UNPK   ZWL16,PWL8
         MVC    CWL16,ZWL16
         OI     CWL16+15,X'F0'
         MVI    CWL16,C'+'
         CP     PWL8,=P'0'
         BNM    *+8
         MVI    CWL16,C'-'
         MVC    CWL80+0(15),CWL16
         MVC    CWL80+9(1),=C'.'   /k  (15-6=9)
         XPRNT  CWL80,80           display a
         L      R13,4(0,R13)
         LM     R14,R12,12(R13)
         XR     R15,R15
         BR     R14
         DS     0F
K        DC     PL8'1000000'       10^6 
A        DS     PL8
G        DS     PL8
AN       DS     PL8
* ****** SQRT   *******************
SQRT     CNOP   0,4                function sqrt(x)
         ZAP    X,PWL8
         ZAP    X0,=P'0'           x0=0
         ZAP    X1,=P'1'           x1=1
WHILE2   EQU    *                  while x0!=x1
         ZAP    PWL8,X0
         SP     PWL8,X1
         CP     PWL8,=P'0'         (x0-x1)!=0
         BE     EWHILE2
         ZAP    X0,X1              x0=x1
         ZAP    PWL16,X
         DP     PWL16,X1
         ZAP    XW,PWL16(8)        xw=x/x1
         ZAP    PWL8,X1
         AP     PWL8,XW
         DP     PWL8,=P'2'
         ZAP    PWL8,PWL8(7)
         ZAP    X2,PWL8            x2=(x1+xw)/2
         ZAP    X1,X2              x1=x2
         B      WHILE2
EWHILE2  EQU    *
         ZAP    PWL8,X1            return x1
         BR     R14
         DS     0F
X        DS     PL8
X0       DS     PL8
X1       DS     PL8
X2       DS     PL8
XW       DS     PL8
* end SQRT
PWL8     DC     PL8'0'
PWL16    DC     PL16'0'
CWL80    DC     CL80' '
CWL16    DS     CL16
ZWL16    DS     ZL16
         LTORG  
         YREGS  
         END    AGM