Newton 法

最近関数プログラミングというのが流行っているわけですが、FORTRAN の父 John Backus は、関数プログラミングに対しても先駆者となっています。
昔 von Neumann モデルの限界がどうのこうの巷で騒がれたころに、Bakus のチューリング賞受賞講演 Can Programming Be Liberated from the von Neumann Style? http://www.stanford.edu/class/cs242/readings/backus.pdf を読んだのですが、まるっきりちんぷんかんぷんでした。
最近の関数プログラミングの流行とともに犬猫バカチョンでも分かると称する解説なども増えてきて、あらためて件の論文を読み返すと、多少はわからなくもないような気がしてきます。
よくネットで引用される文章として John Hughes『なぜ関数プログラミングは重要か』http://www.sampou.org/haskell/article/whyfp.html （原文Why Functional Programming Matters http://www.cse.chalmers.se/~rjmh/Papers/whyfp.pdf)というものがあります。この中の Newton 法の例題を今風の Fortran で解いてみることにします。
近似値列を引数のスタックに積んで再帰呼び出し、収束した時点でその数値列を返り値として出力します。

実行結果

ソース・プログラム

Fortran2008 では内部副プログラムを引数として与えることが許可されます。
なお収束しない場合スタックがあふれて異常終了します。

module m_sub
  implicit none
  
  abstract interface
    pure elemental real function t_f(a, x)
      real, intent(in) :: a, x
    end function t_f
  end interface

contains

  recursive function xiter(f, a, xlist) result(res)
    procedure(t_f) :: f
    real, intent(in) :: a, xlist(:)
    real, allocatable :: res(:)
    real :: f0, f1
    f0 = xlist(size(xlist))
    f1 = f(a, f0)
    if ( abs(f1 - f0) < 0.5 * spacing(f0) ) then
      res = xlist
    else  
      res = xiter(f, a, [xlist, f1])
    end if
    return
  end function xiter
  
end module m_sub

program iteration
  use m_sub
  implicit none
  real :: a
  a = 3.0
  print  '(*(g0))', "Newton's method"
  print '(/*(g0))', 'inverse of ', a 
  print '((7f11.6))', xiter(f, a, [0.1])
  print '(/*(g0))', 'square root of ', a 
  print '((7f11.6))', xiter(g, a, [1.0])
  
  a = 7.0
  print  '(//*(g0))', "Newton's method"
  print '(/*(g0))', 'inverse of ', a 
  print '((7f11.6))', xiter(f, a, [0.1])
  print '(/*(g0))', 'square root of ', a 
  print '((7f11.6))', xiter(g, a, [1.0])
  
  stop
contains
  pure elemental real function f(a, x)  ! 1 / a
    real, intent(in) :: a, x
    f = x * (2 - a * x)
    return
  end function f

  pure elemental real function g(a, x)  ! sqrt(a)
    real, intent(in) :: a, x
    g = 0.5 * (x + a / x)
    return
  end function g
end program iteration