C99 and C++17 both support hexadecimal floating point constants so it should be possible to exactly represent the beeb's 40-bit constants.

BigEd wrote: ↑Fri Apr 27, 2018 3:45 pm

Somehow I'd got the idea that early BBC Basic versions used the ratio of two polynomials, but that seems not to be the case.

You were right. The beeb uses continued fractions in this form:

Nn+x/

...

N3+x/

N2+x/

N1

I tried multiplying this out to see if it I could figure out how to calculate the coefficients. This was a dead end but the sequence you get as you add more terms (working up from the bottom) goes like this:

N2+x/N1

(N1N2N3+(N1+N3)x)/(N1N2+x)

(N1N2N3N4+(N1N4+N1N2+N3N4)x+x^2)/(N1N2N3+(N1+N3)x)

Each extra term flips it over and multiples in another x (and stuff), so it is in fact a ratio of two polynomials. With seven constants N1..N7 the polynomials are both cubic. This is probably a well known fact about continued fractions.

I still wanted to know how to figure out the constants without using Wolfram Alpha.

Euler's method is straightforward but produces continued fractions like this:

c(0)+c(1)x/

1-c(2)x/

1+c(2)x-c(3)x/

1+c(3)x-c(4)x/

1+...

This has variables in both the addends and the numerators. I couldn't find a way to convert this to a beeb-style continued fraction.

So how does Wolfram Alpha produce continued fractions in this form?:

c(0)+c(1)x/

1+c(2)x/

1+c(3)x/

1+...

I found this

olde PhD thesis (PDF). On page 71 algorithm 3.16 produces continued fractions in this form:

c(0)/

1+c(1)x/

1+c(2)x/

1+...

I implemented this [N.B. detail missing from the paper: b(n,0)=a(n) n=1,2,3...] and it works but doesn't directly correspond to the beeb's or Wolfram's version because they expect an undivided x at the top. But this is fixable.

Given a power series:

a(0)+a(1)z+a(2)z^2+a(3)z^3...

use the algorithm from the paper to get a continued fraction for

a(1)+a(2)z+a(3)z^2...

(i.e. discard a(0) and divide by z). Then

a(0)+c(0)x/

1+c(1)x/

1+c(2)x/

1+...

is the Wolfram-style continued fraction, which is easily converted to beeb-style.

I implemented this and it gives exactly the same values as Wolfram Alpha. However, once converted to beeb-style these are slightly different to the values used in BASIC 2. I guess this is because of the Chebychev optimisation but I have no idea how to do that.