Mathematical Curiosities

This post contains some interesting mathematical curiosities I’ve recently learned about.

Formula for Primes

C. P. Willans (1964, paper)

Not a fast formula to compute, but it does provide the n^th prime. It is very inefficient from a computational point of view, with two loops embedded.

Wilson’s theorem: $(j-1)!+1$ is divisible by j precisely when j is prime or 1. Hence, the formula uses it with $cos$ and square to test for primality.

Bertrand’s postulate: for every integer $m \geq 2$ there exists a prime number $p$ such that $m < p < 2m$. Hence the formula sums until $2^n$

However, according to the paper “What is an answer?” by Herbert Wilf, this algorithm is not a real answer to the question “what is the n^th prime?”

Here’s a short Python code implementing the formula.

import math

def prime(n):
    sum = 0
    for i in range(1, int(math.pow(2, n)) + 1):
        inner_sum = 0
        for j in range(1, i + 1):
            div = (math.factorial(j - 1) + 1) / j
            inner_term = math.pow(math.cos(math.pi * div), 2)
            inner_sum = inner_sum + math.floor(inner_term)
        term = math.pow(n / inner_sum, 1 / n)
        sum = sum + math.floor(term)
    return 1 + sum

for i in range(1, 7):
    print(prime(i), end=", ")

The code above correctly prints 2, 3, 5, 7, 11, 13 but fails when i = 7 with a calculation error, and then at i = 8 with a float overflow. I’ve tried using numpy (most notably the numpy.float128 type) and Decimal, but the values reached by the variables in the inner loop are astronomical, and it keeps failing when calculating the next valid values, 17 and 19.

The Ruby code below uses the BigDecimal library and calculates the first 10 values correctly, but it takes around 4 minutes to do so on my laptop:

require 'bigdecimal'
require 'bigdecimal/math'

class BigDecimal
  def fact
    (2..self).reduce(1,:*)
  end
end

PI = BigMath.PI(10)

def prime(n)
  sum = BigDecimal(0)
  limit = 2 ** n
  for i in 1..limit do
    inner_sum = BigDecimal(0)
    for j in 1..i do
      div = (BigDecimal(j - 1).fact + 1) / BigDecimal(j)
      inner_term = BigMath.cos(PI * div, 10) ** BigDecimal(2)
      inner_sum = inner_sum + inner_term.floor
    end
    term = (n / inner_sum) ** Rational(1, n)
    sum = sum + term.floor
  end
  sum + BigDecimal(1)
end

for i in 1..10 do
  puts(prime(i).to_s('F'))
end

These are the times my laptop took to find the n^th prime, and as you can see the times increase dramatically after the 8^th prime! I haven’t tried searching for the 12^th, my estimation is that the calculation would take around 5 hours to complete.

Stirling’s Approximation

This is a jaw-dropping formula, found in the 18th century, involving $\pi$ and $e$ that returns the n^th factorial with increasing accuracy. There’s a page about it on Wikipedia.

Found through:

• Gamma function
• Taylor series
• Gaussian integral

Here’s a table comparing this approximation to the actual values of the factorial function, and the relative error decreases substantially after $n$ reaches 4 or 5:

Download the (source spreadsheet).

Curious Integral

This is the kind of exercises I like, found it on Twitter when it was still an interesting place.

Quake Fast Square Root Algorithm

This is the famous (and wickedly fast) code used in Quake to calculate square roots, which was one of the secrets of the incredible speed of the game.

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y = number;
    i = * ( long * ) &y;                     // evil floating point bit hack
    i = 0x5f3759df - ( i >> 1 );             // what the fuck?
    y = * ( float * ) &i;
    y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    return y;
}

There are blog posts, videos, and a Wikipedia article. I also wrote an article on De Programmatica Ipsum about this algorithm.

Note: the current implementation of sqrt() in math.h is much faster than the one available back in the 90s when Quake was made, so today this code is no longer needed. But back then it was nothing else than a miracle.

Somebody wrote the mandatory version in Rust of the code above, with the same clarification that the current sqrt() implementation is actually faster and more accurate.

Simple Square Root Algorithms in Python

More square root algorithms, but these are definitely not as optimized; these are just using both a recursive and non-recursive version of Newton’s algorithm, for the sake of illustration.

import math

def sqroot1(number, guess = 100, iterations = 20):
    if iterations == 0:
        return guess
    current = (guess + number / guess) / 2
    return sqroot1(number, current, iterations - 1)

def sqroot2(number, precision = .01):
    # Adapted from exercise 2.1 of "Shaum's Computer Science" page 182
    guess = float(number / (1 + len(str(number))))
    iterations = 100
    while abs(guess * guess - number) > precision and iterations > 0:
        iterations -= 1
        guess = (guess + number / guess) / 2
    return guess, iterations

input = 9876543210
val1 = sqroot1(input)
val2, iterations2 = sqroot2(input, 0.00001)
val3 = math.sqrt(input)
print(f'sqroot1({input}) = {val1} (20 iterations)')
print(f'sqroot2({input}) = {val2} ({iterations2} iterations)')
print(f'math.sqrt({input}) = {val3}')

Running the python code above returns this:

sqroot1(9876543210) = 99380.79900061179 (20 iterations)
sqroot2(9876543210) = 99380.79900061179 (82 iterations)
math.sqrt(9876543210) = 99380.79900061178