A blog

blog#A Dumb Way To Calculate Pi.md

Let $f(x) = \sqrt{1 - x ^ 2}$, $f^{(n)}(x)$ be the n-th derivative of $f$.

Based on the Taylor series expansion (around 0), we have:

$$ f(x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(0)}{n!} x^n $$

So,

$$ \int_0^1 f(x) dx = \sum_{n=0}^{\infty} \frac{f^{(n)}(0)}{(n+1)!} $$

Besides,

$$ \begin{aligned} \pi &= 4 \times \int_0^1 \sqrt{1 - x ^ 2} dx \\ &= 4 \times \int_0^1 f(x) dx \\ &= 4 \times \sum_{n=0}^{\infty} \frac{f^{(n)}(0)}{(n+1)!} \end{aligned} $$

The problem now is how to get $f^{(n)}(0)$.

Based on some calculation, we can get:

$$ \begin{aligned} f^{(0)}(x) &= (1 - x ^ 2) ^ \frac{1}{2} \\ f^{(1)}(x) &= - x (1 - x ^ 2) ^ {-\frac{1}{2}} \\ f^{(2)}(x) &= - (1 - x ^ 2) ^ {-\frac{1}{2}} - x ^ 2 (1 - x ^ 2) ^ {-\frac{3}{2}} \\ f^{(3)}(x) &= - 3 x (1 - x ^ 2) ^ {-\frac{3}{2}} - 3 x ^ 3 (1 - x ^ 2) ^ {-\frac{5}{2}} \\ f^{(4)}(x) &= - 3 (1 - x ^ 2) ^ {-\frac{3}{2}} - 18 x ^ 2 (1 - x ^ 2) ^ {-\frac{5}{2}} - 15 x ^ 4 (1 - x ^ 2) ^ {-\frac{7}{2}} \\ \\ \\ ... \end{aligned} $$

Based on the derivatives calculated above, we can observe a pattern in the coefficients. Let's organize them in the following table:

The table shows the coefficients of terms of the form $x^k (1-x^2)^{-\frac{m}{2}}$ in each derivative $f^{(n)}(x)$, where the rows represent the derivative order $n$, and the columns represent the exponent $k$ of $x$ (with $k$ and $m$ satisfying $m = n + k - 1$).

$n$	$k=0$	$k=1$	$k=2$	$k=3$	$k=4$
0	$1$
1		$-1$
2	$-1$		$-1$
3		$-3$		$-3$
4	$-3$		$-18$		$-15$

If you find the derivatives yourself, you may find that any of the coefficient, $k_0$ and $n_0$ for instance, only affects $k_0 - 1$ (if $k_0 - 1 >= 0$) and $k_0 + 1$ in $n_0 + 1$.

For $k_0 - 1$, it adds $coefficient \times k_0$ to it.
For $k_0 + 1$, it adds $coefficient \times (n_0 + k_0 - 1)$ to it.

To get $f^{(n)}(0)$, obviously $f^{(2n + 1)}(x) = 0, n \in \mathbb{N}$,
and $f^{(2n)}(0) n \in \mathbb{N}$ just equals to it's coefficient where $k = 0$.

So we can have the program.

const N: usize = 10000;
const PRECISION: u32 = 256;

use rug::{Float, Integer};

fn main() {
    println!("Calculating pi with {} terms...", N);

    let mut cur_n = 2;
    let mut vec = vec![Integer::new(); N + 10];

    let mut collector: Vec<Integer> = Vec::with_capacity(N / 2 + 10);
    collector.push(Integer::from(1));

    vec[0] = Integer::from(-1);
    vec[2] = Integer::from(-1);

    while cur_n <= N {
        if cur_n % 2 == 0 {
            for i in (1..cur_n + 2).step_by(2) {
                vec[i] = Integer::new();
            }
            for i in (0..cur_n + 1).step_by(2) {
                let k = vec[i].clone();

                vec[i + 1] += (i + cur_n - 1) * &k;
                if i >= 1 {
                    vec[i - 1] += i * k;
                }
            }

            collector.push(std::mem::take(&mut vec[0]));
            // println!("collector added: {}", collector[collector.len() - 1]);
        } else {
            for i in (0..cur_n + 2).step_by(2) {
                vec[i] = Integer::new();
            }
            for i in (1..cur_n + 1).step_by(2) {
                let k = vec[i].clone();

                vec[i + 1] += (i + cur_n - 1) * &k;
                vec[i - 1] += i * k;
            }
        }
        cur_n += 1;
    }

    let mut pi = Float::with_val(PRECISION, 0);
    let mut factorial_base = Integer::from(1);
    let mut factorial_cap = Integer::from(2);
    for i in 0..collector.len() {
        pi += Float::with_val(PRECISION, std::mem::take(&mut collector[i]))
            / Float::with_val(PRECISION, factorial_base.clone());

        factorial_base = factorial_base * &factorial_cap;
        factorial_base = factorial_base * (&factorial_cap + Integer::from(1));
        factorial_cap += 2;
        // println!("pi_{} = {:.16}", i, pi);
    }
    pi *= 4;

    println!("pi = {}", pi);
}

I wrote this in hurry so it's not really optimized. Sorry about that.

Calculating pi with 10000 terms...
pi = 3.141593717260361692706004115729853061820382289134097618662650072056492604360802

That's quite near... I mean, in Astronomy they use 3 as $\pi$...

Just kidding. This infinite series converges slowly, but it uses a quite basic method. I know there're better ways like Trigonometric Functions and Binomial Theorem, but I don't know how to use them. I'm new to calculus.