Re: continued fractions

Jonah Thomas <jethomas5@xxxxxxxxx> wrote:
Hugh Aguilar <hugoaguilar@xxxxxxxxxxxx> wrote:

http://groups.google.co.jp/group/comp.lang.forth/msg/6fa4e33f552ca10a

Note that when you get a continued fraction for pi the issue is not
that you have a continued fraction for pi. What you will do with it is
multiply it by stuff. And the best continued fraction is the one that
gives you the right product the most times, with perhaps some other
properties like you might want the error to be low half the time and
high half the time.

But ignoring all that, if you want to see how good n * a / b is, you
can take the original double-precision x and multiply it by the
single-precision n to get a triple-precision product. Then you divide
by the single-precision b to get a double-precision quotient. Then you
take the low cell and use it to decide whether to round the high cell
up or down. That single-precision result is what you compare n * a / b
by to see whether it's perfect or off by one. If it's off by more than
one then you don't have a very good continued fraction yet.

Check whether I got the details right, it's early for me and I could
have them all wrong. But the basic idea I'm sure is right. After you
do*/ you will have a single-precision result. Is that the best
single-precision result or is it off? To find out, you might as well
start with your high-precision "correct" form and calculate the best
single-precision result that you will compare against. I don't think
this needs DUM/MOD .

Note that */ does not round, it truncates. I think this will throw a
systematic bias into your results. Continued fraction methods give you
the best result with rounding, and this does not round. So something
other than the best continued fraction might be what you want. I could
be wrong about this too.

It seems to me that on average, you could hope for a continued
fraction pair that gives you the best single-precision result about
half the time, and gives you a result that is 1 bit low 1/4 the time
and a result that is 1 bit high 1/4 the time. Or, since there are so
many continued fraction pairs to choose from, maybe you could get one
that gives you the correct result much more than half the time and
still makes the errors high or low about equally. For 16 bits you can
easily do an exhaustive search and find the best pairs completely by
brute force. For 32 bits it's harder....

I don't think I've ever gotten a response on this, but I got interested
in it one more time so I'll write about it again.

-----------
When people create continued fractions the usual approach is to try to
get a pair of numbers x y so that x/y is very close to some target. But
for fixed point math this is a red herring.

For fixed point math, there is a best fixed point approximation. For
example, if you want a 16-bit approximation for 10000*pi then you can't
get any better than 31416. 31415 is clearly not as good and 31417 is
worse. Never mind how close x/y would come to pi if you did an
infinite-precision x/y. When you finish your single-precision arithmetic
you will have a single-precision result and it will be off by however
much it's off. It might be the best possible single-precision
approximation or it might be off by one. If it's off by more than one
that's worse.

So for any x y pair, here's a brute-force method to see how good it is.
I'll talk about pi as an example.

First you find the biggest number you can multiply by pi and still get a
valid result. For signed 16 bits that would be 10430. 10431*pi = 32769+
which is out of bounds.

Then, for every number z from zero to the biggest number, you do x y
*/ and get a value. Then you convert z to floating point, multiply by
pi, round, subtract, and convert back to single-precisition fixed point.
If the result is zero then you have the best possible fixed point
approximation for z*pi. If the result is 1 or -1 then you don't have the
best approximation.

So for every possible z you calculate how far off z x y */ puts you.
Ideally, you would get the best approximation every time. Failing that,
you'd like to have only a few errors and you'd like them to be high
about as often as they're low. The pattern of rounding errors will not
necessarily be smallest when x/y is closest to pi.

-----------

Second point. The Forth */ does not round. Its division truncates. This
is not wrong behavior. The routine would slow down and use more stack
space if it tracked the remainder and adjusted the last bit. If you want
a */ that rounds you can write one without a lot of trouble. I did that
and I never found an example where I really needed it.

However, the concept behind continued fractions and foley fractions
requires rounding. Would you expect to get the same rounding errors when
you truncate instead? There is no particular reason to expect continued
fractions to give you the x y pair that will most often give the best
result with */ .

-----------

I considered rethinking continued fractions to work with truncation. But
instead I tried for a brute force proof of concept. With today's
desktops it's fast to check all of the 16-bit possibilities.

\ -----------------------
\ test the rounding errors

2variable approx \ the current x/y pair
FVARIABLE OK \ the correct value

\ compare single-precision approximation to best approximation
\ where "best" comes from floating point.
: test1 ( n -- -1|0|1)
DUP approx 2@ */ \ n x
SWAP
S>D D>F OK F@ F* FROUND F>D DROP
- ;

: ARRAY
CREATE CELLS ALLOT
DOES> SWAP CELLS + ;

\ I assume error will not be more than 1 bit.
2 ARRAY ERROR
: clear-errors
3 0 do 0 i ERROR ! loop ;

\ test error for every possible multiplier.
\ with symmetric division the negative values will mirror positives.
: tests ( -- )
clear-errors
10431 0 DO
I TEST1 1 SWAP 1+ 0 max 2 min ERROR +!
LOOP ;

: t ( n d -- )
approx 2! tests 0 error ? 1 error ? 2 error ? ;

\ -------------------

355 113 approx 2!
pi OK F!



\ -------------------
\ brute force testing

variable best

\ check whether a particular x y pair is good
\ if it works for one value then test it for all.
: bf
10430 approx 2@ */ 32766 32769 within if
tests 1 ERROR @ best @ > if
1 ERROR @ best !
0 ERROR ? 1 ERROR ? 2 ERROR ?
approx 2@ . . ." "
then
then ;

bfs 5171 5259 0 113 355 1668 8255 507 120 377
1958 8258 214 833 2617 1920 8266 244 953 2994
1892 8270 268 1073 3371 1839 8271 320 1433 4502 ok

355 113 is low almost half the time. I believe this is because */
truncates. By the criteria I tested for, 4502 1433 does better -- it
gives the best possible result almost 80% of the time. But when it gives
the wrong result that result is almost 6 times as likely to be low as
high.

355 113 is wrong almost half the time and the error is highly biased.
4502 1433 is better.



I tried choosing a minimum percentage to be right and then seeing how
well I could balance the errors.

VARIABLE SIDES

: bf'
10430 approx 2@ */ 32766 32769 within if
tests 1 ERROR @ 8000 > if
0 ERROR @ 2 ERROR @ - ABS DUP SIDES @ U< if
SIDES !
0 ERROR ? 1 ERROR ? 2 ERROR ?
approx 2@ . . ." " else
drop
then
then
then ;
\ try all plausible 16-bit x y pairs.
\ sloppy but fast enough.

: bfs'
65535 SIDES !
\ try every x value
32767 4 DO
\ try every plausible y value
I DUP 3 / OVER 4 / ?DO
dup I approx 2! bf'
key? if key 32 = if exit then then
LOOP drop
LOOP ;

bfs' 1668 8255 507 120 377 1425 8029 976 607 1907
1413 8006 1011 1701 5344 1410 8002 1018 2795 8781
1410 8001 1019 9479 29780 ok

I can get less bias on the error, but it results in somewhat bigger
errors.


I instinctively wanted the error to be unbiased. But consider pi.
3141592.... If you did have errors up to half a bit out of 16 bits, that
would go from 3141542 to 3141642, which would round to 31416 92% of the
time and down to 31415 .08% of the time.

If the error was bigger it could go from 3141534 to 3141650 which would
round to 31416 86% of the time and down to 31415 close to 14%. With 8000
correct, we'd have around 7 instances rounding up to 31417 per 23
rounding down to 31415. So maybe that's more the ratio to aim for.

I'm still not completely clear about all this, but I'm clear on some
things. If you're looking for a pair of numbers x y so that z x y */ is
as close as possible to z * OK for real numbers, your criteria should be
about how often you get the best approximation and the bias of the
errors, not about how close z x y */ would come if they were all real
numbers. And */ truncates while continued fractions assume it will
round.
.

Relevant Pages

  • continued fraction and approximation
    ... if I impose the additional constraint a/b<=x? ... Candidates for algorithm: ... Approximate x from below by a continued fraction ... Use normal continued fraction and take last approximation ...
    (sci.math)
  • Re: Status of Waring-problem
    ... that was again a helpful short-course in approximation ... > where ak is the first partial quotient of the continued fraction omitted ... the best we can do with rational approximations is the convergents in the ...
    (sci.math)
  • continued fractions for golden ratio and a coincidence?
    ... I calculated the first few values for the continued fraction ... which turns out to be the continued fraction representation ... we can note that the nth approximation is ... to the continued fraction for the golden ratio, ...
    (sci.math)
  • Re: continued fractions for golden ratio and a coincidence?
    ... which turns out to be the continued fraction representation ... we can note that the nth approximation is ... to the continued fraction for the golden ratio, ... m a positive integer. ...
    (sci.math)
  • Re: continued fractions
    ... 64-bit system to calculate rationals for a 32-bit system. ... you have a continued fraction for pi. ... the single-precision b to get a double-precision quotient. ... Note that */ does not round, ...
    (comp.lang.forth)