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
Wouldn't this be better?
ReplyDeletedefine numdigits(v) {
return length(v)-scale(v)
}
@zsbana: Thanks for the comment on the built-in function length, integrated it.
ReplyDelete