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 a₁, a₂, a₃, … its Dirichlet series generating function (dsgf) is defined to be

The Dirichlet series generating function for the sequence of all 1’s is the Riemann zeta function ζ(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

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

a₁, a₂, a₃, …

and B(s) is the dsgf of

b₁, b₂, b₃, …

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.

Related posts

• Dirichlet convolution
• Computing ζ(3)
• Bernoulli numbers and impossible sums
• Distribution of Riemann zeta zeros

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

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

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

and for n > 1 define

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

then

where

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

• Transforms and convolutions
• Multi-part definitions in LaTeX
• Counting primitive bit strings

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

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

and that the sum is not really infinite because all the terms with index i larger than log₅ 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  f∘f and f∘f∘f 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

• Floor, ceiling, bracket
• Complex floor

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/p^a has? Indeed there is, and usually the answer is

r p^a-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

10^x = 12⁴⁰

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

• Cyclic fractions
• Three quirks of bc
• Last digits of Fibonacci numbers

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

I was reading a book on computer algebra and ran across an interesting theorem about Carmichael numbers in the one of the exercises. I’ll present that theorem below, but first I’ll back up and say what a pseudoprime is and what a Carmichael number is.

Fermat’s theorem

If p is a prime number, then for any integer b, Fermat’s little theorem says

b^p-1 = 1 (mod p).

The contrapositive of this theorem is useful for testing whether a number is (not) prime: if for some b,

b^p-1 ≠ 1 (mod p)

then p is not a prime. For example, we could test whether 1003 is prime by computing [1]

2¹⁰⁰² (mod 1003).

When we do, we find the result is 990, not 1, and so 1003 is not prime.

Pseudoprimes

Let’s test whether 341 is prime. If we compute

2³⁴⁰ (mod 341)

we get 1. But clearly 341 = 11 × 31. If we go back and read the top of the post carefully, we see that Fermat’s theorem only gives a way to prove a number is not prime. We say 341 is a pseudoprime base 2 because it is a composite number that slips past the primality test based on Fermat’s little theorem with b = 2.

Strictly speaking we should say 341 is a Fermat pseudoprime base 2. When someone says “pseudoprime” without further qualification, they usually mean Fermat pseudo prime. But there are other kinds of pseudoprimes. For any theorem that gives a necessary but not sufficient for a number to be prime, the composite numbers that satisfy the condition are called pseudoprimes. See this post for examples, such as Euler pseudoprimes and Wilson pseudoprimes.

Note that 341 is not a Fermat pseudoprime base 3 because

3³⁴⁰ = 56 (mod 341)

If a number passes Fermat’s primality test for several bases, it’s probably prime. (Here’s a blog post that carefully unpacks what “probably prime” means.)

Carmichael numbers

We saw that 341 is a Fermat pseudoprime for base 2 but not base 3; it slipped past our first check on primality but not our second. Some numbers are more slippery. Carmichael numbers are composite numbers that slip past Fermat’s primality test for every base b.

Actually, not every base. A Carmichael number n passes Fermat’s primality test for any base b relatively prime to n.

The smallest Carmichael number is 561. There are infinitely many Carmichael numbers, but they’re distributed very thinly.

Searching for Carmichael numbers

Now I’m finally ready to get to the exercise I mentioned at the top of the post. I’m reading Computer Algebra by Wolfram Koepf. The exercise asks the reader to prove that if for a given k the three numbers 6k+1, 12k+1, and 18k+1 are all prime, then their product is a Carmichael number. This theorem was discovered by J. Chernick in 1939.

We can use this to search for Carmichael numbers with the following Python code. (Note that this approach doesn’t produce all Carmichael numbers—it doesn’t produce 561, for example.) It is unknown whether it produces an infinite sequence of Carmichael numbers.

    from sympy import isprime

    for k in range(2, 100):
        if isprime(6*k+1) and isprime(12*k+1) and isprime(18*k+1):
            print(k, (6*k+1)*(12*k+1)*(18*k+1))

Testing values of k up to 100 reveals six Carmichael numbers.

     6 294409
    35 56052361
    45 118901521
    51 172947529
    55 216821881
    56 228842209

Related posts

• Testing for primes less than a quintillion
• Searching for Fermat primes
• How likely is it that a probable prime is prime?

[1] This may not look like a good idea at first. 2¹⁰⁰² is an enormous number. Surely it would be easier to test whether 1003 is prime than to compute such a large number. Except you don’t have to compute 2¹⁰⁰² per se: you have to compute 2¹⁰⁰² (mod 1003). You can reduce your intermediate results (mod 1003) at every step, so you never have to work with large numbers. You can also use fast exponentiation, and so Fermat’s primality test is much more efficient than it may seem at first glance.

Python’s pow function optionally takes a modulus as a third argument. So, for example,

    pow(2, 1002)

computes 2¹⁰⁰², and

    pow(2, 1002, 1003)

computes 2¹⁰⁰² (mod 1003). It produces the same result as

    2**1002 % 1003

but is more efficient. Not that it matters in this example, but it matters more when the exponent is much larger.

The previous post looked at repeatedly applying the function s(n) which is the sum of the divisors of n less than n. It is an open question whether the sequence

s( s( s( … s(n) … ) ) )

always converges or enters a loop. In fact, it’s an open question of whether the sequence starting with n = 276 converges.

A heuristic argument for why the sequence ought to converge, at least much of the time, is that s(n) is usually less than n. This can be made precise as explained here: in the limit as N goes to infinity, the proportion of numbers less than N for which s(n) < n is at least 0.752.

Even though applying s usually leads to a smaller result, conceivably it could lead to a much larger result when it increases, as with the hailstone problem.

I made some histograms of the ratio s(n) / n to get a feel for how much s increases or decreases its argument. (I imagine this has all been done before, but I haven’t seen it.)

Here’s my first attempt, for n up to 1000.

So apparently the ratio is often less than 1.

I thought it would be clearer to work on a log scale and center the graph at 0 so I next made the following graph.

This shows that the aliquot function s never increases its argument by much, and it often reduces its argument quite a bit, at least for numbers less than 1000.

Next I made a histogram with numbers less than 1,000,000 to see whether the pattern continued. Here’s what I got.

The most obvious feature is the huge spike. If you can ignore that, mentally cropping the image to just the bottom part, it looks basically like the previous graph.

What’s up with that spike? It’s in the bin containing log(0.5). In other words, the ratio s(n)/n is often near 0.5.

So I thought maybe there’s some number theoretic reason the ratio is often exactly 1/2. That’s not the case. It’s only exactly 1/2 once, when n = 2, but it’s often near 1/2.

In any case, applying s very often returns a smaller number than you started with, so it’s plausible that the iterated aliquot sequence would often converge. But empirically it looks like some sequences, such as the one starting at n = 276, diverge. They somehow weave through all the numbers for which the sequence would come tumbling down.

Let s(n) be the sum of the proper divisors of n, the divisors less than n itself. A number n is called excessive, deficient, or perfect depending on whether s(n) – n is positive, negative, or 0 respectively, as I wrote about a few days ago.

The number s(n) is called the aliquot sum of n. The iterated aliquot problem asks whether repeated iterations of s always converge to a perfect number or to 0. This problem has something of the flavor of the Collatz conjecture, a.k.a. the hailstone problem.

As mentioned here, about 1/4 of all integers are excessive and about 3/4 are deficient. That is, s(n) is smaller than n for about 3 out of 4 integers. Based on just this, it seems plausible that iterated aliquots would often converge. (See the next post for a visualization of how much the function s increases or decreases its argument.)

We can investigate the iterated aliquot problem with the following Python code.

    from sympy import divisors
    
    def s(n):
        if n == 0:
            return 0
        return sum(divisors(n)) - n
    
    def perfect(n):
        return n == s(n)
    
    def iterate(n):
        i = 0
        while True:
            if perfect(n) or n == 0:
                return i
            n = s(n)
            i += 1

The function iterate counts how many iterations it takes for repeated applications of s to converge.

For n = 1 through 24, the sequence of iterates converges to 0, the only exception being the perfect number 6.

When n = 25 the sequence converges, but it converges to 6 rather than 0. That is, 25 is the smallest number which is not perfect but for which the sequence terminates in a perfect number.

When n = 138, the sequence converges, but it takes 178 iterations to do so. For smaller n the sequence converges much sooner.

When n = 220, we see something new. s(220) = 284, and s(284) = 220. That is, our code is stuck in an infinite loop.

So we need to go back and refine our question: does the sequence of iterations always converge to 0, converge to a perfect number, or enter a loop?

We have to modify the Python code above to look for loops;

    def iterate(n):
        i = 0
        seen = set()
        while True:
            if perfect(n) or n == 0 or n in seen:
                return i
            seen.add(n)
            n = s(n)
            i += 1

At least for n < 276 the sequence converges.

When I tried 276, the iterations got much larger and slower. I let it run for 400 iterations and the last iteration was

165895515315097460916365008018720.

This took 150 seconds. I switched over to Mathematica and it calculated the same result in 0.6 seconds.

Then I asked Mathematica to run 600 iterations and that took 64 seconds. I went up to 700 iterations and that took 6133 seconds. So it seems that each set of 100 iterations takes about 10 times longer than the previous one.

The last result I got was

19672082490505732759807289983767208175341724164433534581895943000246457396.

Maybe the iterations turn around eventually, but so far it looks like they’ve reached escape velocity. The growth isn’t monotonic. It goes up and down, but more often goes up.

The problem of whether the sequence always converges is unsolved. I don’t know whether it’s known whether the particular sequence starting at 276 converges.

Update: See the first comment. It’s an open problem whether the sequence starting at 276 converges.