Seeking proof for the formula relating Pi with its convergents

Not giving the solution, but some ideas that could lead to it.

If this formula works, it will be because the integral can be decomposed on $P(x) + \frac{1}{1+x^2}$. It's easy to integrate a polynomial (and will give you a rational number), and $1 \over {1+x^2}$ will give you arctan and as such $\pi \over 4$.

So, taking the problem in reverse, you can look at the expansion of arctan between 0 and $\pi \over 4$. You can then convert that polynomial in a continued fraction that will give you the approximations you're looking at.

This reference could help http://www.math.binghamton.edu/dikran/478/Ch7.pdf (p10)

I'll try to have a closer look later.

---

EDIT: Ok, so after peeking through Mathoverflow as well, here is the idea.

Starting from the very general following integral: $I_n = \int_0^1 \frac{x^l(1-x)^m(\alpha + \beta x^2) }{\gamma(1+x^2)}dx$.

As previously discussed, you want to be able to represent this as follows: $\int_0^1 \frac{P(x)(1+x^2) + C}{\gamma(1+x^2)}dx$ with $P(x)$ a polynomial and $C$ a constant.

With some algebra on the polynomials you know that $P(x)= Q(x)(1+x^2) + Ax+B$. So you want to find $A=0$ and if possible get some ideas on $B$.

You need to evaluate the above polynomial on $x=i=e^{i\frac{\pi}{2}}$ and $-i$ to identify the coefficient.

$ A=0 \Leftrightarrow P(i)=P(-i)=0 \Leftrightarrow A = \Im (i^l(1-i)^m(\alpha - \beta) ) = (\alpha - \beta) 2 ^\frac{m}{2} \sin(\frac{l \pi}{2} + \frac{-m\pi}{4} ) = (\alpha - \beta) 2 ^\frac{m}{2} \sin(\frac{\pi}{4} (2l-m) )$

So we have a condition on $2l-m$ to zero $A$ which is $\frac{\pi}{4} (2l-m) = K\pi \Leftrightarrow 2l-m \equiv 0 [4]$.

In particular, $m=2m'$ which we will use going forward and $l-m' \equiv 0[2]$. We can set $l-m' = 2 \epsilon$

Second part is to look at $B$

$B= \Re (i^l(1-i)^m(\alpha - \beta) ) = (\alpha - \beta) 2 ^\frac{m}{2} \cos(\frac{\pi}{4} (2l-m) ) = (\alpha - \beta) 2 ^{m'} \cos(\frac{\pi}{2} (l-m') ) = (\alpha - \beta) 2 ^{m'} (-1)^\epsilon$

So, to sum up, we have $I_n = \int_0^1 Q(x) + \frac{B}{\gamma(1+x^2)}$. As you're trying to approximate $\pi$, you need to take $\gamma = B/4$.

So, provided that $\int_0^1 Q(x) dx$ can be used to approximate fractions of $\pi$, which is likely given the numbers of degress of freedom, we've proved that $I_n$ would be of the following form, which is slightly better as you can drop $j$ from your variables and it shows some relationships better (provided that $\alpha - \beta \not = 0)$:

$(-1)^n (\pi- \frac{p_n}{q_n}) = \int_0^1 \frac{x^{\epsilon+2m'}(1-x)^{2m'}(\alpha + \beta x^2) }{(\alpha - \beta) 2 ^{m'-2} (-1)^{\epsilon}(1+x^2)}dx$.

In your last example, that yields: $n=6,m'=6,\epsilon = -8,\alpha = 77159, \beta = 124360$

I leave it to you to verify it works on the rest of them.

I computed some more values and got convinced this formula isn't relevant yet unless you remove at least the three variables $j$, $l$ and $m$. It is actually very easy to trivially find arbitrarily many tuple $(i,j,k,l,m)$. First of all, if I understand well, your $j$ value seems to be $(m/2-2)$, or am I wrong?

Even after having discarded $j$, it looks like you can find $(i,k)$ solutions for almost any arbitrary tuple $(l,m)$; thus I suggest you find some interesting rule on $(l,m)$ in order to minimize $(i,k)$ if you want your formula to be strong enough. Of course using only one $i$ variable instead of $(i,j,k,l,m)$ would be a significant improvement.

The code below can help you; it is Maxima code yielding thousands and thousand $(i,j,k,l,m)$:

display2d:false$

A002485: [ 3,22,333,355,103993,104348,208341,312689,833719,1146408,4272943,
           5419351,80143857,165707065,245850922,411557987,1068966896,2549491779,
           6167950454,14885392687,21053343141,1783366216531,3587785776203,
           5371151992734,8958937768937 ] $

A002486: [ 1,7,106,113,33102,33215,66317,99532,265381,364913,1360120,
           1725033,25510582,52746197,78256779,131002976,340262731,811528438,
           1963319607,4738167652,6701487259,567663097408,1142027682075,
           1709690779483,2851718461558 ] $

for n:2 thru 26 do (
  for l:1 thru 36 do (
    for m:1 thru 36 do (
      e: expand(integrate(x^l *(1-x)^m*(k+(i+k)*x^2)/(1+x^2), x, 0, 1)),
      myi: coeff(e, i, 1),
      p: coeff(myi, %pi, 1),
      myi2: myi - p*%pi,
      if coeff(myi2, log(2), 1) = 0 then (
        myk: coeff(coeff(e, i, 0),k, 1),
        if myk # 0 then (
          j: round(float(log(abs(p))/log(2))),
          kipos: (-1)^n * (-A002485[n-1]/A002486[n-1])*(2^j) - myi2,
          if kipos # 0 then (
            k1: denom(myk/kipos), i1: num(myk/kipos),
            check1: expand( (k1*myk+p*i1*%pi+myi2*i1)/( abs(i1)*(2^j))*(-1)^(n) - %pi + A002485[n-1]/A002486[n-1] ),
            if check1 = 0 then
              print ("n =",n,"==> (",i1,j,k1,l,m,")" )),
          kineg: -myk/((-1)^n * (-A002485[n-1]/A002486[n-1])*(2^j) + myi2),
          if kineg # 0 then (
            k2: denom(kineg), i2: num(kineg),
            check2: expand( (k2*myk+p*i2*%pi+myi2*i2)/( abs(i2)*(2^j))*(-1)^(n) - %pi + A002485[n-1]/A002486[n-1] ),
            if check2 = 0 then 
              print ("n =",n,"==> (",i2,j,k2,l,m,")" )))))))$

You asked me by mail if I could find some $(i,j,k,l,m)$ solutions for $n=8$, but there is an infinite number of such solutions; for instance:

n = 8 ==> ( -66317 -1 9977 1 2 ) 
n = 8 ==> ( -6963285 3 212651 1 10 ) 
n = 8 ==> ( -66383317 7 833127 1 18 ) 
n = 8 ==> ( -103605391175 11 720252257 1 26 ) 
n = 8 ==> ( -55884747999795 15 517817918873 1 34 ) 
n = 8 ==> ( -66317 0 8805 2 4 ) 
n = 8 ==> ( -2188461 4 89050 2 12 ) 
n = 8 ==> ( -16380299 8 317214 2 20 ) 
n = 8 ==> ( -1305427928805 12 23297755114 2 28 ) 
n = 8 ==> ( 896546759355 14 27371124886 2 32 ) 
n = 8 ==> ( -23756674250925 16 5393931750178 2 36 ) 
n = 8 ==> ( -66317 1 8246 3 6 )
etc.

I took the liberty of using Wolframs Integrator calculator for the following integrals:

(1) $$\int_0^1\frac{x^2(1-x)^2(a+bx^2)}{1+x^2} dx = (\pi-\frac{47}{15})a + (\frac{22}{7} - \pi)b$$ (2) $$\int_0^1\frac{x^2(1-x)^8(a+bx^2)}{1+x^2} dx = (\frac{3959}{315}-4\pi)a + (4\pi - \frac{4838}{385})b$$ (3) $$\int_0^1\frac{x^4(1-x)^4(a+bx^2)}{1+x^2} dx = (\frac{22}{7}-\pi)a + (\pi - \frac{1979}{630})b$$ (4) $$\int_0^1\frac{x^4(1-x)^8(a+bx^2)}{1+x^2} dx = (4\pi - \frac{4838}{385})a + (\frac{566063}{45045} - 4\pi)b$$ (5) $$\int_0^1\frac{x^6(1-x)^4(a+bx^2)}{1+x^2} dx = (\pi - \frac{1979}{630})a + (\frac{10886}{3465} - \pi)b$$ (6) $$\int_0^1\frac{x^6(1-x)^8(a+bx^2)}{1+x^2} dx = (\frac{566063}{45045} - 4\pi)a + (4\pi - \frac{188685}{15015})b$$ (7) $$\int_0^1\frac{x^8(1-x)^4(a+bx^2)}{1+x^2} dx = (\frac{10886}{3465} - \pi)a + (\pi - \frac{141511}{45045})b$$

From this you can see the trend. Since there are 2 degrees of freedom for a,b, with one condition that $a-b = \pm 1$, then their value can be fixed in such a way to match any n-th convergent of $\pi$. This accounts for why there are multiple solutions for $i,j,k,l$ being found. Actually only 2 parameters are necessary.

I don't see the point, since we have only pushed finding the convergent for $\pi$ into finding $a_n,b_n$ for the n-th convergent of $\pi$, any of the above integrals would work, and I am not sure that discovering the rule for generating $a,b$ will be easy.