lambdaway :: agm

7 | list | login |load

| |

_h1 arithmetic-geometric mean

_p From [[rosettacode.org|https://rosettacode.org/wiki/Arithmetic-geometric_mean#Python]] calculate the arithmetic-geometric mean of two numbers.

_h2 python

{pre
def agm(a0, g0, tolerance=1e-10):
    an, gn = (a0 + g0) / 2.0, sqrt(a0 * g0)
    while abs(an - gn) > tolerance:
        an, gn = (an + gn) / 2.0, sqrt(an * gn)
    return an

print agm(1, 1/sqrt(2))  
-> 0.847213084835
}

_p Multi-Precision Version

{pre
from decimal import Decimal, getcontext

def agm(a, g, tolerance=Decimal("1e-65")):
    while True:
        a, g = (a + g) / 2, (a * g).sqrt()
        if abs(a - g) < tolerance:
            return a

getcontext().prec = 70
print agm(Decimal(1), 1 / Decimal(2).sqrt())
-> 0.847213084793979086606499123482191636481445910326942185060579372659734
}

_h2 lambdatalk

{pre
'{def eps 1e-15}
-> {def eps 1e-15}

'{def agm
  {lambda {:a :g}
   {if   {> {abs {- :a :g}} {eps}}
    then {agm {/ {+ :a :g} 2} {sqrt {* :a :g}}}
    else :a }}} 
-> {def agm
  {lambda {:a :g}
   {if {> {abs {- :a :g}} {eps}}
    then {agm {/ {+ :a :g} 2} {sqrt {* :a :g}}}
    else :a }}}

'{agm 1 {/ 1 {sqrt 2}}}
-> {agm 1 {/ 1 {sqrt 2}}}  
}

_p Multi-precision version using the {b lib_BN} library

{require lib_BN}

{pre
'{BN.DEC 70} 
-> {BN.DEC 70}  70 digits
'{def EPS {BN./ 1 {BN.pow 10 45}}}
-> {def EPS {BN./ 1 {BN.pow 10 45}}} 
 
'{def AGM
  {lambda {:a :g}
   {if   {= {BN.compare {BN.abs {BN.- :a :g}} {EPS}} 1}
    then {AGM {BN./ {BN.+ :a :g} 2} {BN.sqrt {BN.* :a :g}}}
    else :a }}}  
-> {def AGM
  {lambda {:a :g}
   {if {= {BN.compare {BN.abs {BN.- :a :g}} {EPS}} 1}
    then {AGM {BN./ {BN.+ :a :g} 2}
              {BN.sqrt {BN.* :a :g}}}
    else :a }}}

'{AGM 1 {BN./ 1 {BN.sqrt 2}}}
-> {AGM 1 {BN./ 1 {BN.sqrt 2}}}
}

_p {i alain marty | 2023/01/11}