Tcl, 100 bytes, 0-indexed

proc P {k n\ 1 d\ 1} {time {set n ([incr d [expr 2*$d*$k]]+$n*$k)
incr k -1} $k
list [expr 2*$n] $d}

Try it online!

Tcl, 102 bytes, 0-indexed

proc P {k n\ 1 d\ 1} {time {set n ([set d [expr 2*$d*$k+$d]]+$n*$k)
incr k -1} $k
list [expr 2*$n] $d}

Try it online!

Tcl, 102 bytes. 1-indexed

proc P {k n\ 1 d\ 1} {while {[incr k -1]} {set n ([incr d [expr 2*$d*$k]]+$n*$k)}
list [expr 2*$n] $d}

Try it online!

Ungolfed

Tcl, 153 bytes, 1-indexed

proc P k {

    set n 1
    set d 1

    while {[incr k -1]} {
        set x [expr {2*$k+1}]
        set d [expr {$d*$x}]
        set n [expr {$d+$n*$k}]
    }

    list [expr 2*$n] $d
}

Try it online!

CASIO BASIC (CASIO fx-9750GIII), 40 bytes

?→N
1→Q
1
Do
Dsz N
2NQ+Q→Q
AnsN+Q
LpWhile N-1
{2Ans,Q

works pretty well

based on l4m2's answer in Node.JS

Haskell, 54 bytes

f n=foldr(\k(a,b)->(k*a+4*k*b+2*b,2*k*b+b))(2,1)[1..n]

0-indexed. Outputs the unreduced fraction as a (numerator, denominator) tuple. Attempt This Online!

Explanation

Uses the same formula as my Zephyr answer, calculated by right-folding over the list [1..n] with a starting value of \$2\$ = (2,1). The fold function is:

$$ \begin{align} g(k, \frac{a}{b}) &= \frac{k}{2k+1} \cdot \frac{a}{b} + 2 \\ &= \frac{ka}{2kb+b} + 2 \\ &= \frac{ka+4kb+2b}{2kb+b}. \end{align} $$

Zephyr, 95 bytes

input k as Integer
set p to 0
while k>0
set p to(p*(k/((2*k)+1)))+2
set k to k-1
repeat
print p

1-indexed. Try it online!

Explanation

Using Zephyr's built-in fractions, we calculate the result in a single variable. By distributing the \$2\$ over the formula, we can get a simpler algorithm:

$$ \begin{aligned} \frac{64}{21} &= 2 \left(1 + \frac{1}{3} \left(1 + \frac{2}{5} \left(1 + \frac{3}{7}\right)\right)\right) \\ &= 2 + \frac{1}{3} \left(2 + \frac{2}{5} \left(2 + \frac{3}{7} \cdot 2\right)\right) \\ &= 2 + \frac{1}{3} \left(2 + \frac{2}{5} \left(2 + \frac{3}{7} \left(2 + \frac{4}{9} \cdot 0\right)\right)\right) \end{aligned} $$

So starting at \$p=0\$, we loop over decreasing values of \$k\$, multiplying \$p\$ by \$\frac{k}{2k+1}\$ and adding \$2\$ at each step.

APL(NARS), 47 chars

r←f w
w+←r←1x
→3×⍳0≥w-←1⋄r←1+r×w÷1+2×w⋄→2
r×←2

It was a little difficult to understand i have to start from the end. But it is very simple... It seems output are all rationals. I failed to write the recursive one.

  ⍪f¨0,⍳15
                        2
                      8r3
                    44r15
                    64r21
                  976r315
               10816r3465
             141088r45045
              47104r15015
           2404096r765765
        45693952r14549535
        45701632r14549535
        80863232r25741485
    5256312832r1673196525
    3153846272r1003917915
457311809536r145568097675
833925152768r265447707525

R, 74 61 bytes

`~`=\(n,m=1)`if`(m>n,1,(a=n~m+1)[1]*(2*m+1)+m*a*0:1)*1:2^!m-1

Attempt This Online!

Recursive approach to construct the Rabinowitz-Wagon formula as stated in the question.

Outputs the (0-based) n-th non-simplified fraction of the sequence as two integers in the order (denominator, numerator).

Ungolfed:

f=function(n){
    g=function(n,m=1){
        if(m>n)c(1,1) else {
            x=g(n,m+1)
            a=x[1];b=x[2]
            c(2*m*b+b+m*a,2*m*b+b)
        }
    }
    g(n)*2:1
}

R, 62 60 54 bytes

\(n){a=1
if(n)for(i in n:1)a=a*i+(T=T*2*i+T)
c(2*a,T)}

Attempt This Online!

Port of shape warrior t's answer

Wolfram Language (Mathematica), 63 45 bytes

Saved 18 bytes thanks to @att

Golfed version. Try it online!

{2,1}Dot@@Array[{{#,a=2#+1},{0,a}}&,#].{0,1}&

Ungolfed version. Try it online!

(*Define the matrix for a given n*)
matrixForN[n_] := {{n, 2   n + 1}, {0, 2   n + 1}}

(*Compute the product of matrices from n=1 to m*)
productOfMatrices[m_] := 
 Fold[Dot, {{2, 0}, {0, 1}}, Table[matrixForN[n], {n, m}]]

(*Apply the final matrix to the vector {0,1}*)
f[m_] := productOfMatrices[m] . {0, 1}

(*Print the results for n=1 to 20*)
Table[{n, f[n]}, {n, 1, 20}] // 
 Do[Print[ele[[1]], " -> ", ele[[2]][[1]], "/", ele[[2]][[2]]], {ele, #}] &

PARI/GP, 28 bytes

n->(z=2)+sum(i=1,n,z/=2+1/i)

Attempt This Online!

A port of @G B's Ruby answer.

\$0\$-indexed. Returns a rational number.

PARI/GP, 43 bytes

n->prod(i=0,n,[i+2*!i,!!i*k=2*i+1;0,k])[,2]

Attempt This Online!

\$1\$-indexed. Returns a column vector of numerator and denominator.

Using the fact that composition of linear fractional transformations (\$x\mapsto\frac{ax+b}{cx+d}\$) corresponds to matrix multiplication.

The \$n\$-th output is \$\begin{bmatrix}2&0\\0&1\end{bmatrix}\begin{bmatrix}1&3\\0&3\end{bmatrix}\cdots\begin{bmatrix}n&2n+1\\0&2n+1\end{bmatrix}\begin{bmatrix}0\\1\end{bmatrix}\$.

R, 84 bytes

\(n,m=0:n*2+1)cbind(Reduce(\(a,b)a*b+2*prod(1:((b-1)/2)),c(2,m),,,T)[-1],cumprod(m))

Attempt This Online!

An anonymous function that takes an integer \$n\$ and outputs the first \$n+1\$ fractions of Rabinowitz-Wagon 𝜋-formula as two columns side by side - the numerator and the denominator, respectively.

## APL, 36 bytes

This is a function that takes n as input and outputs a pair of integers representing the numerator and denominator of the answer:

{{⍺×⍵+⍵[2]0}/(⊂2 1),(⊢,1+2∘×)¨⍳⍵}

APL, 32 bytes

Updated Solution:: ⎕←{⍺×⍵+⍵[2]0}/(⊂2 1),(⊢,1++⍨)¨⍳⎕ - I converted from a function to stdin & stdout, and used suggestion by akamayu (converted 2∘× → +⍨ to save a byte).

Try it online!

I'll try to shorten it by using more tacit style to remove named arguments - any suggestions are welcome 🙂.

Kap: 21 characters

(1+)⍛×/⌽2,÷∘(1+2×)1+⍳

This function takes n, and returns a rational number. This takes advantage of the built-in support for rational numbers in Kap, where dividing integers always yields a rational. To convert the result to floating point, one can add 0.0 to the result.

Ruby, 38 35 bytes

->n{(z=2)+(1..n).sum{|x|z/=2+1r/x}}

Try it online!

M, 11 bytes

RḤ‘İ×R×\S‘Ḥ

Try it online!

Interestingly enough, the Jelly equivalent (which doesn't have rational number support) is only 1 byte longer. The TIO link demonstrates the results for each input \$0\$ to \$n\$.

How it works

RḤ‘İ×R×\S‘Ḥ - Main link. Takes an integer n on the left
R           - Generate the range [1, 2, ..., n]
 Ḥ          - Unhalve; [2, 4, ..., 2n]
  ‘         - Increment; [3, 5, ..., 2n+1]
   İ        - Inverse; [1/3, 1/5, ..., 1/2n+1]
     R      - Range; [1, 2, ..., n]
    ×       - Multiply; [1/3, 2/5, ..., n/2n+1]
       \    - Scan by:
      ×     -   Product; [1/3, 1/3×2/5, ...]
        S   - Sum; 1/3 + 1/3×2/5 + ...
         ‘  - Increment; 1 + 1/3 + 1/3×2/5 + ...
          Ḥ - Unhalve; 2 + 2×1/3 + 2×1/3×2/5 + ...

APL(Dyalog Unicode), 31 2724 21 bytes ^SBCS

It outputs the \$n\$th term with the denominator comes before the numerator. Assumes ⎕io←0. 24 -> 21 bytes thanks to @att.

∊2(×∘(+⍀)⌿⊣@0,⍨¨1+×)⍳

Try it on APLgolf!

Try it on APLgolf! (24bytes)

Try it on APLgolf! (27bytes)

Try it on APLgolf! (31bytes)

Charcoal, 27 bytes

ＮθＩＥ²ΣＥ∨ι⊕θΠ⊕⊞ＯＥθ⎇‹νλν⊗⊕ν¬ι

Try it online! Link is to verbose version of code. Explanation: Directly computes the numerator and denominator, inspired by @Bubbler's last approach.

Ｎθ                          Input as a number
    ²                       Literal integer `2`
   Ｅ                        Map over implicit range
        ι                   Current value
       ∨                    Logical Or
          θ                 Input number
         ⊕                  Incremented
      Ｅ                     Map over implicit range
                θ           Input number
               Ｅ            Map over implicit range
                   ν        Innermost value
                  ‹         Is less than
                    λ       Inner value
                 ⎇          If true then
                     ν      Innermost value
                        ν   Else innermost value
                       ⊕    Incremented
                      ⊗     Doubled
             ⊞Ｏ             Append
                          ι Current (outer) value
                         ¬  Incremented
            ⊕               Vectorised increment
           Π                Take the product
     Σ                      Take the sum
  Ｉ                         Cast to string
                            Implicitly print

Python, 71 70 60 57 bytes

f=lambda n,q=1,p=1:n and f(n-1,q:=q*2*n+q,p*n+q)or(2*p,q)

Attempt This Online!

Outputs the (0-indexed) nth term of the sequence as a (numerator, denominator) tuple.

Ungolfed algorithm:

def f(n):
    numerator, denominator = 1, 1
    # i/(2i + 1) = n/(2n + 1), ..., 3/7, 2/5, 1/3
    for i in range(n, 0, -1):
        # multiply by i/(2i + 1)
        numerator *= i
        denominator *= 2 * i + 1
        # add 1 (p/q -> (p + q)/q = p/q + q/q = p/q + 1)
        numerator += denominator
    return 2 * numerator, denominator

For n=3:

  * 3/7     + 1      * 2/5       + 1       * 1/3        + 1         * 2
1 ----> 3/7 --> 10/7 ----> 20/35 --> 55/35 ----> 55/105 --> 160/105 --> 320/105

Raku (Perl 6) (rakudo), 49 43 bytes

{2+2*[+] [\*] map {.FatRat/(2*$_+1)},1..$_}

Attempt This Online!

Vyxal s, 51 bits^v2, 6.375 bytes

ƛƛd›/;Πd

Try it Online!

Bitstring:

000001100001010101101101100001101101001010101111101

Ports the expanded formula from Bubbler's J answer.

Explained

ƛƛd›/;Πd­⁡​‎‎⁡⁠⁡‏‏​⁡⁠⁡‌⁢​‎‎⁡⁠⁢‏⁠‎⁡⁠⁢⁢‏‏​⁡⁠⁡‌⁣​‎⁠‎⁡⁠⁢⁡‏‏​⁡⁠⁡‌⁤​‎‎⁡⁠⁣‏⁠‎⁡⁠⁤‏‏​⁡⁠⁡‌⁢⁡​‎‎⁡⁠⁢⁣‏‏​⁡⁠⁡‌⁢⁢​‎‎⁡⁠⁢⁤‏‏​⁡⁠⁡‌⁢⁣​‎‏​⁢⁠⁡‌­
ƛ         # ‎⁡Over each n in the range [1, in]:
 ƛ   ;    # ‎⁢  Over each m in the range [1, n]:
    /     # ‎⁣    m divided by
  d›      # ‎⁤    2m + 1
      Π   # ‎⁢⁡  Product of each fraction
       d  # ‎⁢⁢  Doubled
# ‎⁢⁣Summed by the s flag
💎

Created with the help of Luminespire.

JavaScript (Node.js), 43 bytes

f=(n,p=q=1)=>n?f(n-1,p*n+(q*=n-~n)):[p*2,q]

Attempt This Online!

Port of shape warrior t's answer

J, 14 bytes

2#.~1%2+1%]-i.

Attempt This Online!

If we expand the formula entirely, we get 2 + 2*1/3 + 2*1/3*2/5 + 2*1/3*2/5*3/7 + ..., which can be interpreted as a mixed base conversion of [..., 2, 2, 2, 2] in base [..., 3/7, 2/5, 1/3].

J, 19 bytes

2(*>:)/@,1%2+1%1+i.

Attempt This Online!

A straightforward implementation of the formula. Takes an arbitrary-precision integer, constructs the array that looks like 2 1/3 2/5 3/7, and reduces from the right by x * (1 + y).

J, 20 bytes

(>:@+:#.0&=,:2*!)@i.

Attempt This Online!

Tried to be clever with some mixed-base magic:

(1 + 1/3(1 + 2/5)) * 3*5 = 3*5 + 1*5 + 1*2
(1 + 1/3(1 + 2/5(1 + 3/7))) * 3*5*7 = 3*5*7 + 1*5*7 + 1*2*7 + 1*2*3

So the numerator can be computed by evaluating [0!, 1!, 2!, 3!, ...] in mixed base [1, 3, 5, 7, ...]. Then the denominator can also be represented using [1, 0, 0, 0, ...]. Might be useful in languages without rational number support.

K (ngn/k), 24 bytes

2 1*1(|+\|*)/+-\(1-)\!-:

Try it online!

Input n. Return the \$n\$th term as an integer pair.