Rabinowitz-Wagon \$\pi\$ formula

Question

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}$$

Task

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 sequence 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 code-golf, 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

shape warrior t · Accepted Answer · 2024-07-29 04:24:51Z

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

Bubbler · Accepted Answer · 2024-07-29 00:55:25Z

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.

att · Accepted Answer · 2024-07-29 00:49:54Z

6

K (ngn/k), 24 bytes

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

Try it online!

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

answered Jul 29 at 0:49

att

21.3k2 gold badges16 silver badges60 bronze badges

Add a comment |

lyxal · Accepted Answer · 2024-07-29 01:21:17Z

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.

G B · Accepted Answer · 2024-08-01 07:01:40Z

4

Ruby, 38 35 bytes

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

Try it online!

edited Aug 1 at 7:01

answered Jul 29 at 9:45

G B

22.6k1 gold badge20 silver badges53 bronze badges

1

\$\begingroup\$ That's very elegant! \$\endgroup\$
– int 21h Glory to Ukraine
Commented Aug 1 at 19:20

Add a comment |

l4m2 · Accepted Answer · 2024-07-29 01:13:49Z

3

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

edited Jul 29 at 1:13

answered Jul 29 at 1:05

l4m2

26.3k2 gold badges23 silver badges99 bronze badges

Add a comment |

bb94 · Accepted Answer · 2024-07-29 01:43:19Z

3

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

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

Attempt This Online!

edited Jul 29 at 1:43

answered Jul 29 at 0:18

bb94

3,67814 silver badges19 bronze badges

Add a comment |

Neil · Accepted Answer · 2024-07-29 06:38:14Z

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

caird coinheringaahin g · Accepted Answer · 2024-07-29 18:31:51Z

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

alephalpha · Accepted Answer · 2024-08-02 03:08:52Z

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

138 Aspen · Accepted Answer · 2024-08-03 00:38:03Z

3

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

edited Aug 3 at 0:38

answered Aug 2 at 1:28

138 Aspen

6,3211 gold badge4 silver badges13 bronze badges

\$\begingroup\$ -18 \$\endgroup\$
– att
Commented Aug 2 at 20:50

Add a comment |

Dominic van Essen · Accepted Answer · 2024-08-23 10:09:19Z

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

akamayu · Accepted Answer · 2024-07-29 13:39:01Z

2

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)

edited Jul 29 at 13:39

answered Jul 29 at 4:37

akamayu

1,0151 silver badge6 bronze badges

1

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

Add a comment |

int 21h Glory to Ukraine · Accepted Answer · 2024-08-01 19:19:45Z

2

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.

answered Aug 1 at 19:19

int 21h Glory to Ukraine

2,4272 gold badges9 silver badges15 bronze badges

Add a comment |

Elias Mårtenson · Accepted Answer · 2024-08-01 08:48:55Z

1

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.

answered Aug 1 at 8:48

Elias Mårtenson

2391 silver badge4 bronze badges

Add a comment |

Sohang Chopra · Accepted Answer · 2024-08-01 10:47:56Z

1

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

edited Aug 1 at 10:47

answered Aug 1 at 8:19

Sohang Chopra

911 silver badge5 bronze badges

1

\$\begingroup\$ 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 . \$\endgroup\$
– akamayu
Commented Aug 1 at 8:56
\$\begingroup\$ @akamayu Thanks 🙂, edited my answer to use your suggestion \$\endgroup\$
– Sohang Chopra
Commented Aug 1 at 10:48

Add a comment |

Stack Exchange Network

Rabinowitz-Wagon \$\pi\$ formula

Task

Test cases

16 Answers 16

Python, 71 70 60 57 bytes

J, 14 bytes

J, 19 bytes

J, 20 bytes

K (ngn/k), 24 bytes

Vyxal `s`, 51 bits^v2, 6.375 bytes

Explained

Ruby, 38 35 bytes

JavaScript (Node.js), 43 bytes

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

Charcoal, 27 bytes

M, 11 bytes

How it works

PARI/GP, 28 bytes

PARI/GP, 43 bytes

Wolfram Language (Mathematica), 63 45 bytes

R, 74 61 bytes

R, 62 60 54 bytes

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

R, 84 bytes

Kap: 21 characters

APL, 32 bytes

Your Answer

Not the answer you're looking for? Browse other questions tagged
code-golf
math
sequence
rational-numbers
pi
or ask your own question.

Hot Network Questions

Rabinowitz-Wagon \$\pi\$ formula

Task

Test cases

16 Answers 16

Python, 71 70 60 57 bytes

J, 14 bytes

J, 19 bytes

J, 20 bytes

K (ngn/k), 24 bytes

Vyxal s, 51 bitsv2, 6.375 bytes

Explained

Ruby, 38 35 bytes

JavaScript (Node.js), 43 bytes

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

Charcoal, 27 bytes

M, 11 bytes

How it works

PARI/GP, 28 bytes

PARI/GP, 43 bytes

Wolfram Language (Mathematica), 63 45 bytes

R, 74 61 bytes

R, 62 60 54 bytes

APL(Dyalog Unicode), 31 2724 21 bytes SBCS

R, 84 bytes

Kap: 21 characters

APL, 32 bytes

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged code-golfmathsequencerational-numberspi or ask your own question.

Related

Hot Network Questions

Vyxal `s`, 51 bits^v2, 6.375 bytes

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

Not the answer you're looking for? Browse other questions tagged
code-golf
math
sequence
rational-numbers
pi
or ask your own question.