# Rabinowitz-Wagon $\pi$ formula

In 1995, Stanley Rabinowitz and Stan Wagon found an interesting algorithm to generate the digits of $$\\pi\$$ one by one without storing the previous results. The algorithm is called the spigot algorithm. It is based on the following formula:

$$\pi = 2 \left(1 + \frac{1}{3} \left(1 + \frac{2}{5} \left(1 + \frac{3}{7} \left(1 + \cdots \frac{k}{2k+1} \left(1 + \cdots\right)\right)\right)\right)\right)$$

This formula can be derived from the Leibniz formula for $$\\pi\$$ using the Euler's transformation.

(The formula was already well-known before the algorithm was found, but I can't find a name for it. So I call it the Rabinowitz-Wagon $$\\pi\$$ formula as @aeh5040 suggested in a sandbox comment.)

Using this formula, we can approximate $$\\pi\$$ by a sequence of fractions. The first few terms are:

\begin{aligned} 2 &= 2, \\ \frac{8}{3} &= 2 \left(1 + \frac{1}{3}\right), \\ \frac{44}{15} &= 2 \left(1 + \frac{1}{3} \left(1 + \frac{2}{5}\right)\right), \\ \frac{64}{21} &= 2 \left(1 + \frac{1}{3} \left(1 + \frac{2}{5} \left(1 + \frac{3}{7}\right)\right)\right), \\ \dots&\dots \end{aligned}

In this challenge, you need to output this sequence of fractions.

You can output a fraction $$\\frac{a}{b}\$$ in any reasonable format, e.g. a built-in rational type, a pair of integers $$\(a, b)\$$, or a string a/b. You don't need to simplify the fraction. For example, you can output $$\\frac{8}{3}\$$ as $$\\frac{40}{15}\$$.

As with standard challenges, you may choose to:

• Take an integer $$\n\$$ as input and output the $$\n\$$th term of the sequence.
• Take an integer $$\n\$$ as input and output the first $$\n\$$ terms of the sequence.
• Take no input and output the sequence indefinitely.

The indices can be 0-based or 1-based.

This is , so the shortest code in bytes in each language wins.

## Test cases

0 -> 2
1 -> 8/3
2 -> 44/15
3 -> 64/21
4 -> 976/315
5 -> 10816/3465
6 -> 141088/45045
7 -> 47104/15015
8 -> 2404096/765765
9 -> 45693952/14549535
10 -> 45701632/14549535
11 -> 80863232/25741485
12 -> 5256312832/1673196525
13 -> 3153846272/1003917915
14 -> 457311809536/145568097675
15 -> 833925152768/265447707525
16 -> 4725585805312/1504203675975
17 -> 14176771899392/4512611027925
18 -> 524540820979712/166966608033225
19 -> 104908189597696/33393321606645


# Python, 717060 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

New contributor
shape warrior t is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.

# 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.

# Vyxals, 51 bitsv2, 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.

# Ruby, 38 35 bytes

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


Try it online!

• That's very elegant! Commented yesterday

# 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

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

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


Attempt This Online!

# 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


# 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 + ...


# 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}\$$.

# APL(Dyalog Unicode), 312724 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)

• 21: ∊2(×∘(+⍀)⌿⊣@0,⍨¨1+×)⍳
– att
Commented Jul 29 at 13:08

# 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.

# 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, #}] &

• -18
– att
Commented 22 hours ago

## 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.

## 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 🙂.

• Nice answer! 2∘×+⍨ saves a byte. You can add a link to razetime.github.io/APLgolf so others can check the answer more easily. It can count the bytes and generate the post for you . Commented 2 days ago
• @akamayu Thanks 🙂, edited my answer to use your suggestion Commented 2 days ago