{{ Program that approximates pi (Archemedes' constant) using the
{{   Gauss-Legendre, or more specifically, the Salamin-Brent algorithm.
{{   See "http://en.wikipedia.org/wiki/Gauss-Legendre_algorithm"
{{
{{ NOTE: This program computes pi the first 13 digits
{{
{{ Juliet 2004.0
{{ Kevin Albrecht 15 March 2004

$lib

{{ Newton's method for square root
rout sqrt:rat (number:rat)
  def count:rat, last:rat, root:rat

  root <- number

  loop (count < 20) and (((root-last) * (root-last)) > 0.0000000000001)
    last <- root
    root <- (number/root + root) / 2
    count <- count + 1

  return root

rout main
  def A:rat <- 1.0, B:rat <- 1.0 / sqrt (2.0)
  def T:rat <- 0.25, P:rat <- 1.0
  def X:rat, Y:rat, pi:rat

  def i:int <- 10
  
  loop i > 0
    X <- (A+B) / 2.0
    Y <- sqrt (A * B)
    T <- T - P * (A - X) * (A - X)
    A <- X
    B <- Y
    P <- 2 * P

    cond A-B < 0.0000000000001
      i <- 0
      goto EXITLOOP
    cond T <> 0
      pi <- ((A + B) * (A + B)) / (4 * T)
    printn "pi: " & str pi

  #EXITLOOP

  printn "Complete..."