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

## 2 comments:

Wouldn't this be better?

define numdigits(v) {

return length(v)-scale(v)

}

@zsbana: Thanks for the comment on the built-in function length, integrated it.

Post a Comment