An approximation of pi in bc
This blog post shows an approximation of the constant ? (3.14...) implemented as a bc program using Machins formula.
Machins 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 -qHere is the equivalent program in dc:
# 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. */
[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
alternative link download