Continued fraction expansion of PI not accurate?-CodePudding

I'm trying to write a Windows C program (Visual C 2019) to generate the continued fraction expansion of π.

The correct values from WolframAlhpa and OEIS shows:

[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1,

However, my values start deviating after 14

[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3, 3, 2, 1, 3, 2, 1, 19,

Here's a demo of the code running.

#include <stdio.h>
#include <stdint.h>
#include <math.h>

int main()
{
    long double u = 3.1415926535897932384626433832795028841971693993751058209749445923;

    printf("[%lld; ",(unsigned long long)u);

    for (int i = 0; i < 20; i  )
    {
        u = 1.0 / (u - floorl(u));
        printf("%lld, ",(unsigned long long)u);
    }

    return 0;
}

Question

Is the code losing some kind of decimal precision, causing the incorrect values?

CodePudding user response：

Use long double constants. Append an L. Otherwise you code is assigning a double value to u and not a long double one.

// long double u = 3.1415926535897932384626433832795028841971693993751058209749445923;
long double u = 3.1415926535897932384626433832795028841971693993751058209749445923L;

With L, my output was nearly the expected:

[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 1, 1,         My output
[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, ... WolframAlhpa

WolframAlhpa is using precision beyond my long double.

On some implementations (I'm thinking Visual C), double and long double are the same precision in which case appending an L will not change things.

Add printf("%La\n", u); after the u declaration to see what value code is truly using.

long double u = 3.1415926535897932384626433832795028841971693993751058209749445923L;
printf("%La\n", u);

0x1.921fb54442d1846ap 1  // Your output may differ.

It is the continued fraction of that value (e.g. machine pi) code is calculating and not of π.

Aside:

Best to use matching specifiers

//printf("[%lld; ",(unsigned long long)u);
  printf("[%llu; ",(unsigned long long)u);

CodePudding user response：

When you write u = 1.0 / (u - floorl(u));, the first time the loop runs it's going to be

u = 1.0/1415926535897932384626433832795028841971693993751058209749445923 and it will have more digits than current pi. So, for each iteration the number of digits of u keeps increasing and at some point, it exceeds the highest precision possible and the value is going to be rounded off.