In the process of writing the previous post, I ran across the Landau-Ramanujan theorem. It turned out to not be what I needed, but it’s an interesting result.

What portion of the numbers less than N can be written as the sum of the squares of two non-negative integers? Edmund Landau discovered in 1908, and Srinivasa Ramanujan independently rediscovered in 1913, that the proportion is asymptotically

c / (log N)^1/2

where c, the Landau-Ramanujan constant, equals 0.76422….

Let’s see how the number of squares less than 1000 compares to the estimate given by the theorem.

from math import sqrt, log

N = 1000
c = 0.76422
print("Predicted: ", c / sqrt(log(N)))

sumsq = N*[0]
for i in range(N):
    for j in range(N):
        n = i**2 + j**2
        if n < N:
            sumsq[n] = 1
            
print("Exact: ", sum(sumsq)/N)

The predicted proportion is 0.291 and the exact proportion is 0.330.

The script takes O(N²) time to run, so we'd need to do something more clever if we wanted to investigate very large N.

The sum of the first n odd numbers is n². For example,

1 + 3 + 5 + 7 + 9 = 25 = 5²

Here’s another way to say the same thing. Start with the list of positive integers, remove every second integer from the list, and take the cumulative sum. The resulting list is the list of squares.

We begin with

1, 2, 3, 4, 5, …

then have

1, 3, 5, 7, 9, …

and end up with

1, 4, 9, 16, 25, …

Cubes

Now start over with the list of positive integers. We’re going to remove every third integer, take the cumulative sum, then remove every second integer, and take the cumulative sum. The result is a list of cubes.

Here are the steps described above for the positive integers up to 20:

1, 2, 3, 4, 5, 6, 7, 8, 10, 11, …, 20
1, 2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20
1, 3, 7, 12, 19, 27, 37, 48, 61, 75, 91, 108, 127, 147
1, 7, 19, 37, 61, 91, 127
1, 8, 27, 64, 125, 216, 343

General case

In general, given an integer k we remove every kth element, take the cumulative sum. Then we reduce k by 1 and repeat. We keep doing this until k = 1, in which case we stop. The end result will be the list of kth powers.

The result at the top of the post, corresponding to k = 2, has been known since ancient times. The generalization to larger k wasn’t known until Alfred Moessner observed it in 1951. Oskar Perron proved Moessner’s conjecture later the same year. The result has been called Moesnner’s theorem or Moesnner’s magic.

Python implementation

A little Python code will make it easier to play with Moessner’s process for larger numbers.

from numpy import cumsum

def moessner(N, k):
    x = range(1, N+1)
    while k > 1:
        x = [x[n] for n in range(len(x)) if n % k != k - 1]
        # print(x)
        x = cumsum(x)
        # print(x)
        k = k - 1
    return x

To see the intermediate steps, uncomment the print statements.

If we run the following code

for k in [2, 3, 4, 5]:
    print( moessner(25, k) )

we get the following.

[1   4   9  16  25  36  49  64  81 100 121 144 169]
[1   8  27  64 125 216 343 512 729]
[1   16   81  256  625 1296 2401]
[1   32  243 1024 3125]

Related posts

• Easy to guess, hard to prove
• Recamán’s sequence

For a positive real number x, the expression x mod 1 means the fractional part of x:

x mod 1 = x − ⌊ x ⌋.

This post will look at the distributions of nx mod 1 and xⁿ mod 1 for n = 1, 2, 3, ….

The distribution of nx mod 1 is easy to classify. If x is rational then nx will cycle through a finite number of values in the interval [0, 1]. If x is irrational then nx will be uniformly distributed in the interval.

The distribution of xⁿ mod 1 is more subtle. We have three basic facts.

1. If 0 ≤ x < 1, then xⁿ → 0 as n → ∞.
2. If x = 1 then xⁿ = 1 for all n.
3. For almost all x > 1 the sequence xⁿ mod 1 is uniformly distributed. But for some values of x it is not, and there’s no known classification for these exceptional values of x.

The rest of the post will focus on the third fact.

Suppose x = √2. Then for all even n, xⁿ mod 1 = 0. The sequence isn’t uniformly distributed in [0, 1] because half the sequence is piled up at 0.

For another example, let’s look at powers of the golden ratio φ = (1 + √5)/2. For even values of n the sequence φⁿ converges to 1, and for odd values of n it converges to 0. You can find a proof here.

At one time it was not known whether xⁿ mod 1 is uniformly distributed for x = 3 + √2.  (See HAKMEM, item 32.) I don’t know whether this is still unknown. Let’s see what we can tell from a plot.

Apparently sometime around n = 25 the sequence suddenly converges to 0. That looks awfully suspicious. Let’s calculate x²⁵ using bc.

    $ echo '(3 + sqrt(2))^25' | bc -l
    13223107213151867.43024035874391918714

This says x²⁵ mod 1 should be around 0.43, not near 0. The reason we’re seeing all zeros is that all the bits in the significand are being used to store the integer part and none are left over for the fractional part. A standard floating point number has 53 bits of significand and

(3 + √2)²⁵ > 2⁵³.

For more details, see Anatomy of a floating point number.

Now let’s see what happens when we calculate xⁿ mod 1 correctly using bc. Here’s the revised plot.

Is this uniformly distributed? Well, it’s not obviously not uniformly distributed. But we can’t tell by looking at a plot because uniform distribution is an asymptotic property.

We didn’t give the definition of a uniformly distributed sequence until now because an intuitive understanding of the term was sufficient. Now we should be more precise.

Given a set E, a sequence ω, and an integer N, define A(E, ω, N) to be the number of elements in the first N elements of the sequence ω that fall inside E. We say ω is uniformly distributed mod 1 if for every 0 ≤ a ≤ b ≤ 1,

lim_{N → ∞} A([a, b), ω, N) / N = b − a.

That is, the relative port of points that fall in an interval is equal to the length of the interval.

Now let’s go back to the example x = √2. We know that the powers of this value of x are not uniformly distributed mod 1. But we also said that for almost all x > 1 the powers of x are uniformly distributed. That means every tiny little interval around √2 contains a number y such that the powers of y are uniformly distributed mod 1. Now if y is very close to √2 and we plot the first few values of yⁿ mod 1 then half of the values will be near 0. The sequence will not look uniformly distributed, though it is uniformly distributed in the limit.

Years ago I wrote a post Golden powers are nearly integers. The post was picked up by Hacker News and got a lot of traffic. The post was commenting on a line from Terry Tao:

The powers φ, φ², φ³, … of the golden ratio lie unexpectedly close to integers: for instance, φ¹¹ = 199.005… is unusually close to 199.

In the process of writing my recent post on base-φ numbers I came across some equations that explain exactly why golden powers are nearly integers.

Let φ be the golden ratio and ψ = −1/φ. That is, φ and ψ are the larger and smaller roots of

x² − x − 1 = 0.

Then powers of φ reduce to an integer and an integer multiple of φ. This is true for negative powers of φ as well, and so powers of ψ also reduce to an integer and an integer multiple of ψ. And in fact, the integers alluded to are Fibonacci numbers.

φⁿ = F_n φ + F_{n − 1}
ψⁿ = F_n ψ + F_{n − 1}

These equations can be found in TAOCP 1.2.8 exercise 11.

Adding the two equations leads to [1]

φⁿ = F_{n + 1} + F_{n − 1} − ψⁿ

So yes, φⁿ is nearly an integer. In fact, it’s nearly the sum of the (n + 1)st and (n − 1)st Fibonacci numbers. The error in this approximation is −ψⁿ, and so the error decreases exponentially with alternating signs.

Related posts

• Plastic powers
• Golden ratio primes
• Golden integration

[1] φⁿ + ψⁿ = F_n (φ + ψ) + 2 F_{n − 1} = F_n + 2 F_{n − 1} = F_n + F_{n − 1} + F_{n − 1} = F_{n + 1} + F_{n − 1}