New WP34S unit

[Nigel (UK)]

(11-12-2023 11:47 PM)Paul Dale Wrote:  I thought the old algorithm was using the continued fraction expansion and that was having problems.


Pauli

That is correct. The algorithm I'm now using extends the continued fraction approach so that the problems go away! At least, I hope so. Here are the details.

We can write the continued fraction representation for a number \(x\) as
\[ (a_0; a_1, a_2, a_3, \ldots),\]
so that
\[x = a_0+{1\over{a_1+\displaystyle{1\over{a_2+\ldots}}}}.\]
These numbers are easy to generate: take the floor of \(x\), record it, find the reciprocal of the fractional part, and repeat. The fractions corresponding to each stage of the continued fraction can be generated recursively from this list of numbers: if the fraction is \(h/k\), then
\[h_n = a_nh_{n-1}+h_{n-2},\quad k_n = a_nk_{n-1}+k_{n-2}.\]
For example, the representation for \(\pi\) is \((3; 7, 15, 1, 292,\ldots)\), leading to fractions
\[3, {22\over 7}, {333\over 106}, {355\over 113},\ldots\]
These are called convergents: each has the property that it is closer to \(\pi\) than any other fraction with the same or smaller denominator.

I believe that the existing WP34S algorithm uses this method. When DENMAX=100 it returns 22/7. However, this isn't the fraction closest to \(\pi\) with denominator less than 100; that best fraction is 311/99. So the WP34S algorithm is not perfect, although it is short, fast, and never wildly wrong.

The C47 algorithm improves on this by bracketing the number \(x\) by a pair of rational numbers, \(a/b\) and \(c/d\). It repeatedly calculates the mediant \((a+c)/(b+d)\) of these numbers, updates the bracket, and repeats the process. This works well for "awkward" numbers like \(\pi\), but is really slow for numbers like 0.5 or 0.3333333333. This lack of speed is an issue on the WP34S.

What I've done is to use the idea of "semiconvergents" discussed in this Wikipedia page. Suppose you reach a stage where the denominator of a convergent is greater than DENMAX. Instead of going back to the previous convergent you can instead generate "semiconvergents" which have denominators intermediate in value. Each of these semiconvergents is a best rational approximation to \(x\), in the sense that it is closer to \(x\) than any fraction with the same or smaller denominator. My algorithm chooses the semiconvergent with the largest denominator less than or equal to DENMAX: in the case of \(\pi\) and DENMAX=100, this is 311/99.

A semiconvergent is the fraction
\[{h\over k}={mh_{n-1}+h_{n-2}\over mk_{n-1}+k_{n-2}},\]
where \(m\) is an integer between \(a_n/2\) and \(a_n\). It's easy to find the value of \(m\) that gives the largest denominator less than DENMAX, so finding the best semiconvergent is fast. There are a few conditions on \(m\), described on the Wikipedia page, but the whole process is straightforward and seems to work well.

Wikipedia doesn't include proofs of the assertion that these semiconvergents are best rational approximations, but there is such a proof here. I've scanned it but I haven't worked through it in detail.

I like this method because, assuming that Wikipedia is correct, it gives the best rational approximation in a relatively small number of steps. In particular it doesn't crash the WP34S hardware, which is a big point in its favour!

Nigel (UK)