I'm having trouble finding a loop invariant to prove one of the algorithms from Knuth's The Art of Computer Programming, Volume 1. Can anyone help me with this?
Here's he algorithm:
Assume L[k] is a table of log_b (2^k/(2^k -1)) for k from 1 to whatever...
Here are the rules for calculating y = log_b(x) for 1 <= x < 2
L1. y <- 0; z <- x/2; k <- 1
L2. if x=1, stop (answer is y)
L3. if x-z < 1, goto L5
L4. x <- x-z; z <- x / 2^k; y <- y + L[k]; goto L2
L5. z <- z/2; k <- k + 1; goto L2
where A <- B means B is assigned to A
ie, in C++ for base 2 logs it could be:
Code:
namespace math {
const double L2[32]= { 1.0000000000000000000,
0.41503749927884381855,
0.19264507794239589256,
0.093109404391481470676,
0.045803689613124791194,
0.022720076500083529651,
0.011315313227834146720,
0.0056465631411420624219,
0.0028205190623786626940,
0.0014095702546713535408,
0.00070461297658937274157,
0.00035226347162902138342,
0.00017612098427402405727,
0.000088057804580026382116,
0.000044028230441777209006,
0.000022013947263955502017,
0.000011006931643385186350,
5.5034553246245392772598580437741e-6,
2.75172503805526697221048749368255e-6,
1.37586186296463416424364914705656e-6,
6.8793076746672366935775758198988e-7,
3.4396534272948303367605503515722e-7,
1.7198266111377426060931611323475e-7,
8.599132799414562175013392308985e-8,
4.29956633563874719242616679057e-8,
2.149783151802240599790735457643e-8,
1.074891571896837110458251460864e-8,
5.37445784947347765328405357126e-9,
2.6872289222340618612134241283e-9,
1.3436144604913616904149611786093319e-9,
6.718072300892635353052178399137866e-10,
3.359036150055274401952526040858497e-10
};
template <class T> T abs (T x) {
return x<0? -x:x;
}
const double DOUBLE_TOL=0.00000000000725;
double log2 (double x)
{
double y = 0.0;
double z = x/2.0;
int k=1;
double pow_k = 2.0;
while (math::abs(x-1.0000) >= math::DOUBLE_TOL) {
double t=x-z;
if (t<1.00000) {
z /= 2.0;
++k;
pow_k *= 2.0;
} else {
x = t;
z = x/pow_k;
y += math::L2[k-1];
}
}
return y;
}
};
Anyone have a good idea how to prove its correctness? It's not for a class or anything, I'm just stumped trying to find a loop invariant.
