Turning the Golay problem sideways

I’ve written a couple posts about the Golay problem recently, first here then here. The problem is to find all values of N and n such that

S(N, n) = \sum_{k=0}^n \binom{N}{k}

is a power of 2, say 2p. Most solutions apparently fall into three categories:

  • n = 0 or n = N,
  • N is odd and n = (N-1)/2, or
  • n = 1 and N is one less than a power of two.

The only two other solutions as far as I know are the two that Golay found, namely N = 23 with n = 3 and N = 90 with n = 2. Those are the only solutions for N < 30,000.

Solutions for fixed p

Phil Connor wrote me with a good suggestion: for fixed n, S(N, n) is an nth degree polynomial in N, and we can search for when that polynomial equals 2p.

Instead of looping over N and seeing whether any S(N, n) is a power of 2 comes out, we could fix p and see whether there’s a corresponding value of N.

Now

S(N, n) = (1/n!)Nn + … + 1

and so we are looking for positive integer solutions to

Nn + … + n! = n! 2p.

Using the rational root theorem

By the rational root theorem, a necessary condition for

Nn + … – n!(2p -1)= 0

to have a rational root is that N must divide the constant term

-n!(2p -1).

That means that for fixed n and fixed p there are only a finite number of Ns that could possibly be solutions, namely the divisors of

n!(2p -1).

Say we wanted to find all solutions with p = 32 and n = 3. Then we only need to check 96 values of N because 3!(232 -1) has 96 divisors as the following Python code shows.

    >>> from sympy import divisors
    >>> len(divisors(6*(2**32-1)))
    96

We could evaluate S(N, 3) for each of the divisors of 6(232 -1) and see if any of them equal 232. And since S(N, n) is an increasing function of N, we could speed things up by being a little clever in our search. For example, if S(d, 3) > 232 for a particular divisor d, there’s no need to keep trying larger divisors. More generally, we could do a binary search.

Not only could this approach be more efficient, I suspect it casts the problem in a form that would make it easier to prove theorems about.

Using the quadratic equation

The rational root theorem isn’t the only theorem we could apply to the polynomial formulation of the Golay problem. For example, if n = 2 we can use the quadratic equation and look at the discriminant to show that a necessary condition for there to be an integer solution is that

2p+3 – 7

must be a positive integer. When p = 12 we get

215 – 7 = 1812,

and this corresponds to Golay’s solution N = 90. There cannot be any solutions for larger values of p if p is odd because then 2p+3 is a perfect square and there cannot be another perfect square so close. I don’t know whether there can be larger even values of p with

2p+3 – 7

being a perfect square. If someone could show there are no such solutions, this would prove Golay found the only solution with n = 2. Maybe something similar would work for larger values of n.

Rumors of a proof

The conjecture I’ve been writing about is listed as Theorem 9.3 in Martin Erickson’s book Introduction to Combinatorics, Erickson does not give a proof but refers to Vera Pless’ book Introduction to the Theory of Error-correcting Codes, Wiley, 1989. However, people on MathOverflow who have seen Pless’ book say there’s no proof there.

Terry Tao says he doubts there is a proof:

To be honest, I doubt that this “theorem” is within the reach of existing technology. … My guess is that Erickson misread the discussion in Pless’s book.

Part of the confusion is that exceptional solutions to the problem discussed here are necessary (but not sufficient) for the existence of perfect Golay codes. It is widely reported that there is only one perfect binary Golay code, one based on S(23, 3) = 211. Some have taken that to mean that there are no more exceptional solutions, but this doesn’t follow. Maybe there are more exceptional solutions, but they don’t lead to the existence of more perfect binary Golay codes.

I haven’t seen a proof that Golay’s 23-bit code is the only perfect binary Golay code, though I’ve seen this theorem quoted repeatedly. It would be good to track down the proof of this theorem (if there is one!) and see whether it proves what I’m calling the Golay conjecture as an intermediate step. (Update: it does not. Seems Terry Tao was right to suspect there is no proof.)

Binomial sums and powers of 2

Marcel Golay noticed that

\binom{23}{0} + \binom{23}{1} + \binom{23}{2} + \binom{23}{3} = 2^{11}

and realized that this suggested it might be possible to create a perfect code of length 23. I wrote a Twitter thread on this that explains how this relates to coding theory (and to sporadic simple groups).

This made me wonder how common it is for cumulative sums of binomial coefficients to equal a power of 2.

Define

S(N, n) = \sum_{k=0}^n \binom{N}{k}

Golay’s observation could be expressed as saying

S(23, 3) = 211.

How often is S(N, n) a power of 2?

There are two easy cases. First,

S(N, 0) = 20

for all N because there’s only one way to not choose anything from a set of N items. Second,

S(N, N) = 2N

for all N. (Proof: apply the binomial theorem to (1 + 1)N.)

If N is odd, then we have

S(N, (N-1)/2) = 2N-1

by symmetry of the binomial coefficients.

Also,

S(2p – 1, 1) = 2p

In summary, S(N, k) is a power of 2 if

  • k = 0 or k = N,
  • N is odd and k = (N-1)/2, or
  • k = 1 and N is one less than a power of two.

Are there any more possibilities? Apparently not many.

The only exceptions I’ve found are (23, 3) and (90, 2). Those are the only possibilities for N < 10,000. (Update: Searched up to 30,000.)

Golay found both of these solutions and mentioned them in [1]. In Golay’s words, “A limited search has revealed two such cases.” He didn’t say how far he looked, and apparently he didn’t have a proof that these were the only exceptional solutions.

(This post has been edited several times since I first posted it, to reflect corrections and a better understanding of the problem.)

Related posts

[1] Golay, M. (1949). Notes on Digital Coding. Proc. IRE. 37: 657.

The congruent number problem

A positive integer n is said to be congruent if there exists a right triangle with area n such that the length of all three sides is rational.

You could always choose one leg to have length n and the other to have length 2. Such a triangle would have area n and two rational sides, but in general the hypotenuse would not be rational.

The smallest congruent number is 5. You can verify that 5 is congruent because the triangle with legs 20/3 and 3/2 has hypotenuse 41/6 and area 5. You can find a list of congruent numbers in OEIS A003273.

You can show that a number is congruent by demonstrating a rational right triangle with the specified area, as with 5 above. Showing that a number is not congruent is harder. How do you know you haven’t looked long enough?

There’s no simple classification of congruent numbers, but there are partial results. Jerrold Tunnell proved [1] a way to show that a number is not congruent.

It suffices to look at square-free numbers, numbers not divisible by the square of an integer. This is because if n = mk² then n is congruent if and only if m is.

Tunnell gives two necessary conditions for whether a square-free number n is congruent: one for the case of n odd and one for the case of n even.

Because Tunnell’s theorem gives necessary conditions, it can only prove that a number is not congruent.

It would be tidy if Tunnell’s conditions were also sufficient, i.e. that they could prove that a number is congruent. And they may be sufficient. Tunnell proved that if the Birch and Swinnerton-Dyer conjecture is true, then the converse of his theorem is true. If you could prove that the Birch and Swinnerton-Dyer conjecture, you could prove the converse of Tunnell’s theorem as a corollary. You could also win a million dollar prize.

Now for Tunnell’s conditions. If an odd square-free positive integer n is congruent, then there are twice as many integer solutions to

2x² + y² + 8z² = n

as

2x² + y² + 32z² = n.

And if n is even and square free, there are twice as many solutions to

4x² + y² + 8z² = n

as

4x² + y² + 32z² = n.

In general it’s hard to know how many integer solutions an equation has, but for the equations above we can find the solutions by brute force because all the terms are positive. We know each variable is between -√n and √n. We could even be more efficient and check smaller ranges of x and z because of the coefficients. And we can mostly look for solutions in the first quadrants and multiply by 8, but we need to be careful not to count solutions with zero components more than once.

Let’s prove that 11 is not a congruent number.

How many solutions to

2x² + y² + 8z² = 11

are there?

Clearly z can only possibly be ±1 or 0. If z = ±1 then x = ±1 and y = ±1. If z = 0 then the only solutions are x = ±1 and y = ±3. So we have 14 solutions: eight of the form (±1, ±1, ±1) and two of the form (±1, ±3, 0).

For the equation

2x² + y² + 32z² = 11

z must be 0, and so we have the eight solutions we found before: all combinations of (±1, ±3, 0).

Our first equation has 14 solutions and the second has 4, and so 11 must not be congruent.

 

[1] Tunnell, Jerrold B. (1983), “A classical Diophantine problem and modular forms of weight 3/2”, Inventiones Mathematicae, 72 (2): 323–334, doi:10.1007/BF01389327

 

Estimating the number of integer partitions

Last week I wrote a few posts that included code that iterated over all partitions of a positive integer n. See here, here, and here. How long would it take these programs to run for large n?

For this post, I’ll focus just on how many partitions there are. It’s interesting to think about how you would generate the partitions, but that’s a topic for another time.

Ramanujan discovered an asymptotic formula for p(n), the number of partitions of n. As n increases,

p(n) \sim \frac{1}{4n\sqrt{3}} \exp\left( \pi \sqrt{\frac{2n}{3}} \right)

The ratio of the two sides of the equation goes to 1 as n → ∞, but how accurate is it for the sizes of n that you might want to iterate over in a program?

Before answering that, we need to decide what range of n we might use. Let’s be generous and say we might want to look at up to a trillion (1012) partitions. How big of a value of n does that correspond to? The estimate above shows p(200) is on the order of 4 × 1012, so we only need to be concerned with n ≤ 200. (See the next post for how to exactly solve for a n that gives a specified value.)

Here’s a plot to show the relative error in the approximation above.

Here’s the Mathematica code that made the plot.

    approx[n_] := Exp[Pi Sqrt[2 n/3]]/(4 n Sqrt[3])
    DiscretePlot[ 1 - approx[n]/PartitionsP[n], {n, 200}]

So the estimate above is always an under-estimate, at least about 4% low over the range we care about. It’s good enough for a quick calculation. It tells us, for example, to be very patient if we’re going to run any program that needs to iterate over all partitions of an integer even as small as 100.

The function PartitionsP above returns the number of partitions in its argument. It returned quickly, so it certainly did not generate all partitions and count them. Algorithms for computing p(n) might make an interesting blog post another time.

Even though it would be impractical to iterate over the partitions of larger values, we still might be curious what the plot above looks like for larger arguments. Let’s see what the plot would look like for n between 200 and 2000.

The shape of the curve is practically the same.

The relative error dropped by about a factor or 3 when we increased n by a factor of 10. So we’d expect if we increase n by another factor of 10, the relative error would drop about another factor of 3, and indeed that’s what happens: when n = 20,000, the relative error is about -0.003.

Partitions and Pentagons

This post will present a remarkable theorem of Euler which makes a connection between integer partitions and pentagons.

Partitions

A partition of an integer n is a way of writing n as a sum of positive integers. For example, there are seven unordered partitions of 5:

  • 5
  • 4 + 1
  • 3 + 2
  • 3 + 1 + 1
  • 2 + 2 + 1
  • 2 + 1 + 1 + 1
  • 1 + 1 + 1 + 1 + 1

Note that order does not matter. For instance, 4 + 1 and 1 + 4 count as the same partition. You can get a list of partitions of n in Mathematica with the function IntegerPartitions[n]. If we enter IntegerPartitions[5] we get

     {{5}, {4, 1}, {3, 2}, {3, 1, 1}, {2, 2, 1}, 
      {2, 1, 1, 1}, {1, 1, 1, 1, 1}}

We say a partition is distinct if all the numbers in the sum are distinct. For example, these partitions of 5 are distinct

  • 5
  • 4 + 1
  • 3 + 2

but these are not

  • 3 + 1 + 1
  • 2 + 2 + 1
  • 2 + 1 + 1 + 1
  • 1 + 1 + 1 + 1 + 1

because each has a repeated number.

We can divide distinct partitions into those having an even number of terms and those having an odd number of terms. In the example above, there is one odd partition, just 5, and two even partitions: 4 + 1 and 3 + 2.

Pentagonal numbers

The pentagonal numbers are defined by counting the number of dots in a graph of concentric pentagons: 1, 5, 12, 22, …. This is OEIS sequence A000326.

The jth pentagonal number is

Pj = j(3j – 1) / 2.

We can define negative pentagonal numbers by plugging negative values of j into the equation above, though these numbers don’t have the same geometric interpretation. (I think they do have a geometric interpretation, but I forget what it is.)

Euler’s pentagonal number theorem

Leonard Euler discovered that the number of even distinct partitions of n equals the number of odd distinct partitions, unless n is a pentagonal number (including negative indices). If n is the jth pentagonal number, then the difference between the number of even and odd distinct partitions of n equals (-1)j. In symbols,

D_e(n) - D_o(n) = \left\{ \begin{array}{ll} (-1)^j & \mbox{if } n = P_j \\ 0 & \mbox{otherwise} \end{array} \right.

Here De is the number of even distinct partitions and Do is the number of odd distinct partitions.

Euler discovered this by playing around with power series and noticing that the series involved were generating functions for pentagonal numbers. Jordan Bell has written a commentary on Euler’s work and posted it on arXiv.

Mathematica examples

Let’s look at the partitions of 12, the 3rd pentagonal number.

IntegerPartions[12] shows there are 77 partitions of 12. We can pick out just the partitions into distinct numbers by selecting partition sets that have no repeated terms, i.e. that have as many elements as their union.

    Select[IntegerPartitions[12], Length[Union[#]] == Length[#] &]

This returns a list of 15 partitions:

    {{12},      {11, 1},   {10, 2},      {9, 3},    {9, 2, 1}, 
     {8, 4},    {8, 3, 1}, {7, 5},       {7, 4, 1}, {7, 3, 2}, 
     {6, 5, 1}, {6, 4, 2}, {6, 3, 2, 1}, {5, 4, 3}, {5, 4, 2, 1}}

Of these sets, 7 have an even number of elements and 8 have an odd number of elements. This matches what Euler said we would get because

7 – 8 = (-1)3.

We’d like to eplore this further without having to count partitions by hand, so let’s define a helper function.

    f[list_] := Module[
        {n = Length[Union[list]]}, 
        If [n == Length[list], (-1)^n, 0]
    ]

This function returns 0 for lists that have repeated elements. For sets with no repeated elements, it returns 1 for sets with an even number of elements and -1 for sets with an odd number of elements. If we apply this function to the list of partitions of n and take the sum, we get

De(n) – Do(n).

Now

    Total[Map[f, IntegerPartitions[12]]]

returns -1. And we can verify that it returns (-1)j for the jth pentagonal number. The code

    Table[Total[Map[f, IntegerPartitions[j (3 j - 1)/2]]], {j, -3, 3}]

returns

    {-1, 1, -1, 1, -1, 1, -1}

We can compute De(n) – Do(n) for n = 1, 2, 3, …, 12 with

    Table[{n, Total[Map[f, IntegerPartitions[n]]]}, {n, 1, 12}]

and get

    {{1, -1}, {2, -1}, {3,  0}, {4,  0}, 
     {5,  1}, {6,  0}, {7,  1}, {8,  0}, 
     {9,  0}, {10, 0}, {11, 0}, {12, -1}}

which shows that we get 0 unless n is a pentagonal number.

Euler characteristic

Euler’s pentagonal number theorem reminds me of Euler characteristic. Maybe there’s a historical connection there, or maybe someone could formalize a connection even if Euler didn’t make the connection himself. In its most familiar form, the Euler characteristic of a polyhedron is

VE + F

where V is the number of vertices, E is the number of edges, and F is the number of faces. Notice that this is an alternating sum counting o-dimensional things (vertices), 1-dimensional things (edges), and 2-dimensional things (faces). This is extended to higher dimensions, defining the Euler characteristic to be the alternating sum of dimensions of homology groups.

We can write De(n) – Do(n) to make it look more like the Euler characteristic, using Dk(n) to denote the number of distinct partitions of size k.

D0(n) – D1(n) + D2(n) – …

I suppose there’s some way to formulate Euler’s pentagonal number theorem in terms of homology.

Dirichlet series generating functions

A couple days ago I wrote about Dirichlet convolution, and in that post I said

Define the function ζ to be the constant function 1.

This was following the notation from the book I quoted in the post.

Someone might question the use of ζ because it is commonly used to denote the Riemann ζ function. And they’d be onto something: Calling the sequence of all 1’s ζ is sort of a pun on the ζ function.

Dirichlet series generating functions

Given a sequence a1, a2, a3, … its Dirichlet series generating function (dsgf) is defined to be

A(s) = \sum_{n=1}^\infty \frac{a_n}{n^s}

The Dirichlet series generating function for the sequence of all 1’s is the Riemann zeta function ζ(s).

\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}

The definition of ζ as a sequence was alluding to its associated Dirichlet series.

Dirichlet convolution

Recall from the earlier post that the Dirichlet convolution of two sequences is defined by

(f*g)(n) = \sum_{d|n} f(d) \,\,g\left(\frac{n}{d} \right)

A basic theorem connecting Dirichlet series and Dirichlet convolution says that if A(s) is the dsgf of

a1, a2, a3, …

and B(s) is the dsgf of

b1, b2, b3, …

then the dsgf of a*b, the Dirichlet convolution of the two series, is A(s) B(s).

In short, the dsgf of a convolution is the product of the dsgfs. This is an example of the general pattern that transforms turn convolutions into products. Different kinds of transformations have different kinds of convolution. See examples here.

Dirichlet series of an inverse

In the earlier post we defined the inverse of a sequence (with respect to Dirichlet convolution) to be the sequence you convolve it with to get the sequence δ(n) given by

1, 0, 0, 0, …

The dsgf of δ(n) is simply the constant 1. Therefore if a*b = δ, then AB = 1. Said another way, if a is the inverse sequence of b, then the dsgf of a is the reciprocal of the dsgf of b.

The earlier post defined the Möbius function μ(n) to be the inverse of the sequence of 1’s, and so its dsgf is the reciprocal of the zeta function.

\sum_{n=1}^\infty \frac{\mu(n)}{n^s} = \frac{1}{\zeta(s)}

Related posts

Dirichlet convolution

That can’t be right

I was skimming a book on number theory [1] recently and came across the following theorem that makes no sense out of context:

An arithmetical function f has an inverse if and only if f(1) ≠ 0.

Wait, what?! That’s obviously false unless I’ve missed something. Maybe “arithmetical function” is far more restrictive than I thought.

No, in this book an arithmetical function is simply a function on positive integers. In general it can be complex-valued, though early in the book all the examples are integer or rational valued.

What I was missing was the meaning of “inverse.” My first thought was that it meant inverse with respect to composition, but it meant inverse with respect to convolution.

Dirichlet convolution

Given to arithmetic functions f and g, their Dirichlet convolution f*g is defined by

(f*g)(n) = \sum_{d|n} f(d) \,\,g\left(\frac{n}{d} \right)

The sum is over all divisors of n and so, for example, the value of f * g at 6 would be

f(1)g(6) + f(2)g(3) + f(3)g(2) + f(6)g(1)

It’s clear that convolution is commutative, and with a little more effort you can show that it’s associative.

There is an identity element for convolution, the function δ(n) defined to be 1 if n = 1 and 0 otherwise. For any arithmetical function f, clearly f*δ = f because all the terms in the sum defining (f*δ)(n) are zero except the term f(n).

Inverse

In the context of the theorem above, the inverse of a function f is another function g whose convolution with f gives the identity δ, i.e. g is the inverse of f if f*g = δ. Inverses, if they exist, are unique.

Now that we have the right definition of inverse in mind, the theorem above isn’t obviously false, but it’s not clear why it should be true. And it’s a remarkable statement: with just the slightest restriction, i.e. f(1) ≠ 0, every arithmetical function has an inverse with respect to convolution. Not only that, the proof is constructive, so it shows how to compute the inverse.

For an arithmetical function f with f(1) ≠ 0, define its inverse function g by

g(1) =\frac{1}{f(1)}

and for n > 1 define

g(n) = -\frac{1}{f(1)} \sum_{\stackrel{d > 1}{d|n}} f(d)\,\,g\left(\frac{n}{d} \right )

This definition is not circular even though g appears on both sides. On the right side g is only evaluated at arguments less than n since the sum restricts d to be greater than 1. The next post implements this recursive definition of g in Python.

Möbius inversion formula

Define the function ζ to be the constant function 1. Since ζ(1) is not zero, ζ has an inverse. Call that inverse μ.

If f = g*ζ, then convolution of both sides with μ shows that f*μ = g. This is the Möbius inversion formula.

When presented this way, the Möbius inversion formula looks trivial. That’s because all the hard work has been moved into prerequisites. Stated from scratch, the theorem is more impressive. Without using any of the definitions above, Möbius inversion formula says that if f is defined by

f(n) = \sum_{d|n}g(d)

then

g(n) = \sum_{d|n}f(d) \, \mu\left(\frac{n}{d} \right)

where

\mu(n) = \left\{ \begin{aligned} {1} &\mbox{\quad if } {n=1} \\ {0} &\mbox{\quad if } {a^2 \,|\, n \mbox{ for some } a > 1} \\ {(-1)^r} &\mbox{\quad if } {n \mbox{ has } r \mbox{ distinct prime factors}} \end{aligned} \right.

We initially defined μ implicitly as the inverse of the constant function 1. When written out explicitly we have the definition of the Möbius function μ above.

Related posts

[1] An Introduction of Arithmetical Functions by Paul J. McCarthy, Springer-Verlag, 1986.

Floor exercises

The previous post explained why the number of trailing zeros in n! is

\sum_{i=1}^\infty \left\lfloor \frac{n}{5^i} \right\rfloor

and that the sum is not really infinite because all the terms with index i larger than log5 n are zero. Here ⌊x⌋ is the floor of x, the largest integer no greater than x.

The post gave the example of n = 8675309 and computed that the number of trailing zeros is 2168823 using the Python code

    sum(floor(8675309/5**i) for i in range(1, 10))

Now suppose you wanted to calculate this by hand. You would want to be more clever than the code above. Instead of dividing n by 5, and then by 25, and then by 125 etc. you’d save time by first computing

⌊8675309/5⌋ = 1735061,

then computing

⌊1735061/5⌋ = 347012,

and so forth.

Is this calculation justified? We implicitly assumed that, for example,

n/25⌋ = ⌊⌊n/5⌋ /5⌋.

This seems reasonable, so reasonable that we might not think to check it, but calculations with floor and ceiling functions can be tricky.

There’s a theorem from Concrete Mathematics that can help us.

Let f(x) be any continuous, monotonically increasing function with the property that whenever f(x) is an integer, x is an integer. Then

f(x) ⌋ =⌊ f(⌊x⌋) ⌋

and

f(x) ⌉ =⌈ f(⌈x⌉) ⌉.

Note that the requirements on f compose. That is, if a function f satisfies the hypothesis of the theorem, so do the compositions  ff and fff etc. This means we can apply the theorem above iteratively to conclude

f(f(x)) ⌋ =⌊ f (⌊f(⌊x⌋)⌋) ⌋

and

f( f(x)) ⌉ =⌈ f (⌈f(⌈x⌉)⌉) ⌉

and similarly for higher compositions.

The hand calculation above is justified by applying the iterated theorem to the function f(x) = x/5.

Related posts

Factorial random number generator

Here’s a strange pseudorandom number generator I ran across recently in [1].

Starting with a positive integer n, create a sequence of bits as follows.

  1. Compute n! as a base 10 number.
  2. Cut off the trailing zeros.
  3. Replace digits 0 through 4 with 0, and the rest with 1.

You’d want to use a fairly large n, but let’s use n = 20 for a small example. Then

20! = 2432902008176640000.

Lopping off the trailing zeros gives

243290200817664

Replacing small digits with 0 and large digits with 1 gives

000010000101110

This is not a practical PRNG, but it does produce a sequence of bits that has good statistical properties.

Several questions naturally arise that will be addressed below.

How many digits does n! have?

I give two ways to compute the answer in this post.

As an example, if we let n = 8675309, then n! has 56424131 digits.

How many trailing zeros does n! have?

An equivalent question is what is the largest power of 10 that divide n!. When we write out the integers from 1 to n, there are more multiples of 2 than multiples of 5, so powers of 5 are the limiting resource when looking to divide n! by powers of 10.

Every fifth number is a multiple of 5, so the number of powers of 5 dividing n! is at least the floor of n / 5. But every 25th number contributes an extra factor of 5. And every 125th contributes still another factor, etc. So the number of factors of 5 in n!, and thus the number of factors of 10 in n!, is

\sum_{i=1}^\infty \left\lfloor \frac{n}{5^i} \right\rfloor

This is formally an infinite sum, but the sum is actually finite because all the terms are zero once the index i is so large that

5i > n.

Or in other words, the number of terms in the sum is bounded by log5 n. Now

log 8675309 = 9.9…

and so the “infinite” sum above has only nine non-zero terms when n = 8675309. A quick Python calculation

    sum(floor(8675309/5**i) for i in range(1, 10))

shows that n! has 2168823 trailing zeros. So n! has roughly 56 million digits, and we chop off roughly the last 2 million digits, giving us 54 million digits, which leads to 54 million bits.

Leading digits

The leading digits are a little biased, because the leading digits of factorials follow Benford’s law. When we carry out the procedure at the top of this post, the first digit of n! will be 1, 2, 3, or 4 more often than it will be 5, 6, 7, 8, or 9. The first bit will be biased, the second bit less biased, the third bit even less biased, etc. The bias declines rapidly, and so would have little effect when generating millions of bits.

In practice you could discard the first few bits, though in practice you’d never use this generator in the first place because it’s so inefficient.

***

[1] Random Numbers and Computers by Ronand Kneusel, Springer, 2018.

Reciprocals of prime powers

Let p be a prime number. This post explores a relationship between the number of digits in the reciprocal of p and in the reciprocal of powers of p.

By the number of digits in a fraction we mean the period of the decimal representation as a repeating sequence. So, for example, we say there are 6 digits in 1/7 because

1/7 = 0.142857 142857 …

We will assume our prime p is not 2 or 5 so that 1/p is a repeating decimal.

If 1/p has r digits, is there a way to say how many digits 1/pa has? Indeed there is, and usually the answer is

r pa-1.

So, for example, we would expect 1/7² to have 6×7 digits, and 1/7³ to have 6×7² digits, which is the case.

As another example, consider

1/11 = 0.09 09 09 …

Since 1/11 has 2 digits, we’d expect 1/121 to have 22 digits, and it does.

You may be worried about the word “usually” above. When does the theorem not hold? For primes p less than 1000, the only exceptions are p = 3 and p = 487. In general, how do you know whether a given prime satisfies the theorem? I don’t know. I just ran across this, and my source [1] doesn’t cite any references. I haven’t thought about it much, but I suspect you could get to a proof starting from the theorem given here.

What if we’re not working in base 10? We’ll do a quick example in duodecimal using bc.

    $ bc -lq
    obase=12
    1/5
    .2497249724972497249

Here we fire up the Unix calculator bc and tell it to set the output base to 12. In base 12, the representation of 1/5 repeats after 4 figures: 0.2497 2497 ….

We expect 1/5² to repeat after 4×5 = 20 places, so let’s set the scale to 40 and see if that’s the case.

    scale = 40
    1/25
    .05915343A0B62A68781B05915343A0B62A6878

OK, it looks like it repeats, but we didn’t get 40 figures, only 38. Let’s try setting the scale larger so we can make sure the full cycle of figures repeats.

    scale=44
    1/25
    .05915343A0B62A68781B05915343A0B62A68781B0

That gives us 40 figures, and indeed the first 20 repeat in the second 20. But why did we have to set the scale to 44 to get 40 figures?

Because the scale sets the precision in base 10. Setting the scale to 40 does give us 40 decimal places, but fewer duodecimal figures. If we solve the equation

10x = 1240

we get x = 43.16… and so we round x up to 44. That tells us 44 decimal places will give us 40 duodecimal places.

Related posts

[1] Recreations in the Theory of Numbers by Alfred H. Beiler.