An approximation of pi in bc

This blog post shows an approximation of the constant π (3.14...) implemented as a bc program using Machin's formula.

Machin's formula is documented on http://en.wikipedia.org/wiki/Numerical_approximations_of_%CF%80#Machin-like_formulae. The implementation below computes acot using its Taylor series expansion, always rounding down for divisions. The implementation knows the number of good digits by computing an upper bound of the total rounding error.

#! /usr/bin/bc -q
# by pts@fazekas.hu at Sat Jul 17 17:38:25 CEST 2010

/** Return an approximation and lower bound for acot(x) * 10 ^ u. */
define acot(x, u) {
  auto sum, xpower, xx, n, term
  sum = xpower = u / x
  xx = x * x
  n = 3
  for (;;) {
    xpower /= xx
    term = xpower / n
    if (term == 0)
      return sum
    sum -= term
    n += 2
    xpower /= xx
    term = xpower / n
    if (term == 0)
      return sum
    sum += term
    n += 2
  }
}

/* Return an integer upper bound for the >= 0 error value pi * 10 ^ u - f(u),
 * where u is 10 ^ nd, f(u) is the integer 4 * (4 * acot(5, u) - acot(239, u)).
 *
 * The magic constants in the maxerr implementation were derived from analyzing
 * the acot implementation, taking into account the rounding (truncation) done
 * in each division.
 */
define maxerr(nd) {
  return (286135312 * nd + 41739380) / 10000000
}

/* Return a string of at most nd characters, prefix of pi,
 * assuming nd >= 4.
 */
define pi(nd) {
  auto u, y
  u = 10 ^ nd
  y = 4 * (4 * acot(5, u) - acot(239, u))
  y /= 10 ^ (length(maxerr(nd)) - 1)
  while (y % 10 == 0) {
    y /= 10
  }
  return y / 10
}

/* Print pi with increasing precision infinitely (until aborted). */
define pinfinite() {
  auto b, nb, nd, a, na
  print "3."
  b = 3
  nb = 1
  nd = 8
  for (;;) {
    a = pi(nd)
    na = length(a)
    /* Print digits not printed yet in previous iterations. */
    print a % (10 ^ (na - nb))
    b = a
    nb = na
    nd *= 3
  }
}

pinfinite()  /* Infinite loop. */

Here is the equivalent program in dc:

[lulx/dspsslxd*sw3snlDx]sC[ls2Q]sS[lplw/dspln/dst0=Slslt-ssln2+snlplw/dsp
ln/dst0=Slslt+ssln2+snlDx]sD[ly10/dsy10%0=Q]sQ[10ld^sz5sxlzsulCx4*239sxlz
sulCx-4*10ld286135312*41739380+10000000/Z1-^/dsy10%0=Qly10/]sP[10lksdlPxd
sgZdsilj-^lgr%nlgshlisjlk3*sklIx]sI[3.]n3sh1sj8sklIx