2

I am using 2021b. I wanted to look at the convergence rate of a iteration algorithm for Pi, using continued fraction. The formula was for wikipedia https://en.wikipedia.org/wiki/Pi#Continued_fractions.

by iteration 10 steps，precision up to 6 digits after dot

 17steps，      10 digits

 18steps，      13 digits

 19steps，      500+ digits

the same above 20 steps.

It seems like Matlab saw through my intention and showed me the "exact" value directly...

This should be a bug, since I tried the same algorithm in both python and scheme, the behaviors are normal.

Now the question is, why is this happenning? I never tried to touch the "exact" value of pi in the code. (the code is below, very short)

    % calculate pi using continued fraction formula
    % pi=4/(1+1^2/(3+2^2/(5+3^2/7+...
    % e.g. take N=15, pi=3.141592653606706

    function p=calc_pi(N)
        p=conti_frac(0,N);
        p=vpa(p,500);
    end
  
    function p=conti_frac(k,N)
        if k==0
            p=(0+4/conti_frac(1,N));
        elseif k==N
            p=2*k-1;
        else
            p=2*k-1+k*k/conti_frac(k+1,N);
        end
    end
Brendan Mitchell
  • 29
  • 5
N.Li
  • 47
  • 7
  • 1
    The question is, why are you using `vpa` there? What is your intention? If you want to show mode than 5 digits of `p`, you can use `format long`, or explicitly call `fprintf` with the right formatting string. – Cris Luengo Aug 11 '23 at 13:08
  • This is embarrassing. I used this coz this is the first thing come to my mind. Theoretically, this should not work, I just tried, and got an incredible result... Thank you. – N.Li Aug 11 '23 at 13:25
  • Oh, I remember what I was thinking. I was expecting my p be a fraction instead of a floating point, under which case what I did made sense. – N.Li Aug 11 '23 at 13:57
  • If you want that, you’d have to use symbolic math inside the recursive function. `vpa` would then evaluate that expression to your chosen number of digits. – Cris Luengo Aug 11 '23 at 14:04
  • Yes, I realized that the moment you told me the key was at vpa. – N.Li Aug 11 '23 at 16:16

1 Answers1

4

You’re doing computations using doubles (64-bit floating point), so you will never be able to compute more than 15 or 16 digits. There just aren’t any more in a double.

What probably happens is that the vpa function decides that the input is so close to pi, you must mean pi, and it substitutes the actual value of pi.

You can verify this is the case by writing

vpa(3.14159265358979, 500)

Add and remove digits to find out how many digits you need for the function to recognize the value of pi.

The documentation to vpa describes this behavior.

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
  • Thank you. I don't know this before. 3.14159265358979 is the exact one matlab wants. – N.Li Aug 11 '23 at 13:28
  • Related discussion here (see comments on last post): https://www.mathworks.com/matlabcentral/answers/549723-why-euler-s-number-e-2-71828-is-not-a-built-in-constant-in-matlab?s_tid=srchtitle – James Tursa Aug 11 '23 at 16:48
  • This is not correct. You will never get more than 15.9 *significant integer digits.* What happens in the fractional part is a completely different thing. For example, 2^-50 can be represented exactly, and there are a lot more than 15 decimal digits in that fraction. – user207421 Aug 14 '23 at 04:16
  • @user207421 I don’t see how why you say contradicts anything I said. – Cris Luengo Aug 14 '23 at 04:47
  • It contradicts 'you will never be able to compute more than 15 or 16 digits'. Try it. – user207421 Aug 14 '23 at 05:26
  • @user207421 Being able to represent any power of two exactly doesn’t help to compute more than 15 or 16 digits of pi, or any other number for that matter. – Cris Luengo Aug 14 '23 at 06:16
  • I suggest you have a look at the decimal representation of 2^-50. It is a counterexample to your claim that 'you will *never* be able to compute more than 15 or 16 digits'. The point is that the available fractional precision depends on the actual value. It isn't *a priori* limited to 15-16 decimal digits. – user207421 Aug 15 '23 at 05:44
  • @user207421 Sure, I could have said "there aren't more than 53 binary digits in a double, which usually corresponds to 15 to 16 decimal digits, but if you're really lucky your number can be represented exactly." But that's a silly thing to say, especially since OP is computing pi, which cannot be represented exactly in any base. In general, you cannot assume you have more than 15 or 16 decimal digits, unless you know your number is one of the special cases. And if you do know that, you don't need this answer. – Cris Luengo Aug 15 '23 at 14:44
  • What you *said* was not 'in genera' but 'never', and I have provided a counterexample which proves your claim false. It's as simple as that. – user207421 Aug 15 '23 at 22:44