problems in accuracy assessment

Dear friends I tried continued fractiona approximation to pi, and got the following rational approximation

1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158

but if I try to find the difference of it with pi, I failed, since

>> sym(1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158)

I got pi directly. How can I assess the difference? Thanks

[John D'Errico]

Please, learn to use the code button to allow your code to be displayed properly. Otherwise, it is unreadable.

[Dingyu Xue]

Thanks John, but it seems I cannot use CODE

Dear friends I tried continued fractiona approximation to pi, and got the following rational approximation

1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158

but if I try to find the difference of it with pi, I failed, since

sym(1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158)

I got pi directly. How can I assess the difference?

Thanks

Answer 1

[John D'Errico]

Using my HPF toolbox, it looks to be correct to about 120 significant digits.

double(log10(abs(hpf('1244969988745789040106366155256015114976454553337906033890313',1000)/hpf('396286255419907262612262286579004636343912658949984351074158',1000) - hpf('pi',1000))))
ans =
    -119.53

Looking back at what you did, this is wrong. Why?

sym(1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158)

ans =
pi

It is wrong because MATLAB will first compute the ratio as a DOUBLE precision number, and only then pass it to sym. sym then recognizes the ratio as a double precision approximation to pi, so it returns that as the result.

If you insist on doing this with the symbolic toolbox, we might try...

log10(abs(vpa(sym('1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158 - pi'),200)))
ans =
-119.52655990941722562288816970381

Of course, if you want a better approximation, say accurate to 250 digits, HPF can provide it too.

[N,D] = rat(hpf('pi',500),hpf('1e-250'))

N =
165213461336854892880967397129207248425612593828416835132438559586109735260365863184040108685528782183687397179966550074377225
D =
52589078074164381557662496030246690228630477463437042437123678996605078304650603819026850043803035007564529614801553827009746

If you need convincing that this ratio is indeed accurate to 250 digits for pi, try this:

double(log10(abs(N/D - hpf('pi',1000))))
ans =
      -250.05

[Dingyu Xue]

Thanks John

Excellent!

Actually I did got the rational number with continued fraction, contfrac in MUPAD, with 120 levels. I tried your method without HPF Toolbox, and made continued fraction with 500 levels, the accuracy is about 500 digits. Is there any hidden rule in the continued fractions? The number of levels is almost the same as the number of correct digits?

[John D'Errico]

Well, you cannot have more correct digits in a rational fraction approximation to a number than are in the original number. So if you are working with 120 digit approximations to pi, you can get roughly 120 correct digits. In the HPF tool, you can provide a tolerance for the approximation. If none is provided, it will assume a tolerance.

And no, in terms of rules, if you are asking if there is any rule for the approximation itself, I don't think that there can exist any simple rule for a rational fraction approximation to a number like pi, although I do not know if that has ever been proven. I'm sure Google would tell me though. Here is the continued fraction I generate for a 100 digit approximation to pi.

3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15 + 1/(-3 + 1/(2 + 1/(2 + 1/(2 + 1/(2 + 1/(3 + 1/(-85 + 1/(-3 + 1/(2 + 1/(15 + 1/(3 + 1/(14 + 1/(-5 + 1/(-2 + 1/(-6 + 1/(-6 + 1/(-100 + 1/(3 + 1/(2 + 1/(6 + 1/(3 + 1/(6 + 1/(-2 + 1/(-6 + 1/(-9 + 1/(9 + 1/(-3 + 1/(-3 + 1/(-8 + 1/(4 + 1/(-2 + 1/(-13 + 1/(3 + 1/(-5 + 1/(2 + 1/(9 + 1/(-2 + 1/(-3 + 1/(8 + 1/(-2 + 1/(-5 + 1/(-2 + 1/(-2 + 1/(-4 + 1/(3 + 1/(4 + 1/(4 + 1/(17 + 1/(-162 + 1/(-46 + 1/(24 + 1/(-3 + 1/(-3 + 1/(6 + 1/(-3 + 1/(-25 + 1/(4 + 1/(-5 + 1/(4 + 1/(-2)))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))