[openbeos] Re: The Scientific Method aka Bresenham NewsletterArticle

  • From: Michael Noisternig <michael.noisternig@xxxxxxxxx>
  • To: openbeos@xxxxxxxxxxxxx
  • Date: Fri, 13 Jun 2003 22:07:15 +0200

Michael Noisternig wrote:

Actually I have not yet been able to find an estimation for the maximum error of those subsequent additions.

Well, it was not that hard to find afterall.
The maximum value we get by all subsequent additions is abs(x2-x1). This value takes m1 = floor(log2(abs(x2-x1)))+1 bits. That leaves m-m1 for the digits after the comma, where m is the number of digits of the mantisse.
If we only use m-m1 bits for the slope "a" from first place then there won't be any additional rounding errors and we get a total error of <= abs(x2-x1)*e', where "e'" is <= 2^-(m-m1+1) for round-to-nearest. This total error is certainly larger than if we hadn't restricted "a" to m-m1 bits. Since 2^m1 < abs(x2-x1)*2 we get
a'-a < abs(x2-x1)^2*e ,
where e is 2^-m.


This is very interesting as we can see that the error for the *addition variant* grows *quadratic* as opposed to the *linear function variant* with a *linearly* growing error.

Solving the bound equation as in the previous example gives us through
        abs(x2-x1)^3*e < 1/2
the answer
        abs(x2-x1) < 2^21
for floating points with 64-bit mantisse.

That means with *some tweaking* the addition code *works* for *integers < 2^21*!

For the sake of continuity some pseudo-code example again (for the 1. octant):

a' = (y2-y1)/(x2-x1);
if ( a' < a )
  a' += 2^-64;  // see next paragraph why
y = 0;
for (x = x1; x < x2; x++) {
  pixel(x,y1+round(y));
  y += a';
}

For checking wether a' </> a one has to think about how to do this as x*a' gets rounded and will likely yield x*a again. One way would be to always subtract/add the maximum error to a' so that a' will always be smaller/greater than a. This will double the actual maximum error however and thus reduce the maximum coordinate distance by 1/sqrt(2).

There is simple answer: the status word of the FPU tells you whether the result was up- or downrounded.



Some pseudo-code example again (for the 1. octant):


a' = (y2-y1)/(x2-x1) + 2^-63;
for (x = 0; x < x2-x1; x++)
  pixel(x1+x*a',y+y1);

Of course that line should have meant pixel(x1+x,y1+round(x*a'));

I wonder if anyone reads my stuff???

Michael Noisternig



Other related posts: