This site is supported by donations to The OEIS Foundation.

User:Peter Luschny/VonMangoldtTransformation

From OeisWiki
Jump to: navigation, search

von Mangoldt

This is a sequel of this blog.

KEYWORDS: Dirichlet operation, Dirichlet transformation, Möbius transformation, von Mangoldt function, von Mangoldt transformation, Chebyshev psi function, least common multiple.

Concerned with sequences: A003418, A056040, A180000, A205957, A205959, A216152, A216153.

Preliminary remark. Implementing number-theoretic functions in the style R = f(A), where A is a list can easily lead to subtle errors. The reason is: The indexing of lists is different in different languages. For instance in Sage (Python) the lists are 0-based, in Maple they are 1-based.

   A = [1,2]; A[0];
   Maple: Error, invalid subscript selector.
   Sage : 1

However, number-theoretic functions are by definition 1-based, so to say. To avoid the danger of 'off-by-one' errors the best strategy is not to pass in lists but functions. Therefore we will switch in this section the signature of the transformations to this format.

The return value R will be still a (Python-)list and as such 0-based. However, if you access R[0] the value returned will be 'undef'. This corresponds to the fact that in the OEIS database the resulting sequences have offset 1.

The Dirichlet operation

The cornerstone of our exposition here is an operation called after Johann Peter Gustav Lejeune Dirichlet.

Dirichlet product

# param:  F, G functions
# param:  n > 0 length
# return: F * G
def dirichlet_product(F, G, n) :
    R = ['undef']
    for k in (1..n) :
        R.append(add(F(d)*G(k/d) for d in divisors(k)))
    return R

Dirichlet prime product

# param:  F, G functions
# param:  n > 0 length
# return: F * G
def prime_divisors(n) : 
    return filter(is_prime, divisors(n)) 

def dirichlet_primeproduct(F, G, n) :
    R = ['undef']
    for k in (1..n) :
        R.append(add(F(d)*G(k/d) for d in prime_divisors(k)))
    return R    

Dirichlet transformation

# param:  F function
# param:  n > 0 length  
# return: transform of [F(1)..F(n)] 
def dirichlet_trans(F, n) :
    R = ['undef']
    for k in (1..n) :
        R.append(add(F(d) for d in divisors(k)))
    return R

Dirichlet prime transformation

# param:  F function
# param:  n > 0 length
# return: transform of [F(1)..F(n)]
def prime_divisors(n) : 
    return filter(is_prime, divisors(n)) 
def dirichlet_primetrans(F, n) :
    R = ['undef']
    for k in (1..n) :
        R.append(add(F(d) for d in prime_divisors(k)))
    return R

Dirichlet 2nd prime transformation

# param:  F function
# param:  n > 0 length
# return: transform of [F(1)..F(n)]
def dirichlet_primetrans2(F, n) :
    R = ['undef']
    for k in (1..n) :
        R.append(add(F(k/d) for d in prime_divisors(k)))
    return R           

In the examples below the following functions are used:

def iota(n) : return 1
def o(n)    : return 0
def nu(n)   : return n 
def oi(n)   : return n % 2
def io(n)   : return 1 - n % 2
def eps(n)  : return 0^(n-1)
def mu(n)   : return moebius(n)
F dirichlet
iota A000005A001221A001221
io A183063A000035 A106404
oi A001227A005087A205745
nu A000203A008472A069359
eps A000012A000004A010051
mu A063524 -A001221A143519

Moebius transformation

# param:  F function
# return: transform of [F(1)..F(n)]
def moebius_trans(F, n) :
    return dirichlet_product(F, mu, n) 

Moebius prime transformation

# param:  F function
# return: transform of [F(1)..F(n)]
def moebius_primetrans(F, n) :
    return dirichlet_primeproduct(F, mu, n) 
F moebius
iota(n) = 1 A063524 A143519
nu(n) = n  A000010 A137851

In our setup the inverse of the Möbius transformation is the Dirichlet transformation and the inverse of the Möbius prime transformation is the Dirichlet prime transformation (already defined in the last section above).

Von Mangoldt's function goes prime

An important example of the Dirichlet product is dirichlet_product(mu, log). Here is what we get:

def lambda(n) :
    return add(moebius(d)*log(n/d) for d in divisors(n))
[round(lambda(n),4) for n in (1..20)]

0.0, 0.6931, 1.0986, 0.6931, 1.6094, 0.0, 1.9459, 0.6931, 1.0986,
0.0, 2.3979, 0.0, 2.5649, 0.0, 0.0, 0.6931, 2.8332, 0.0, 2.9444

Uhh, these are no integers! But wait a moment. If we apply the exponential to the last sequence we get:

1, 2, 3, 2, 5, 1, 7, 2, 3, 1, 11, 1, 13, 1, 1, 2, 17, 1, 19, 1

So now we can easily identify that we are talking about the von Mangoldt function in it's OEISed form A014963. An immediate observation are the many 1s appearing in this sequence. What do they signify? Looking at the indices we see

1, 6, 10, 12, 14, 15, 18, 20, 21, 22, 24, 26, 28, 30, 33, 34, ..

These are the numbers that are not powers of primes, A024619 prepended by 1. Well, this was only the warming up. Next we look at the prime case of the same Dirichlet product (with sign changed for convenience).

def lambda_prime(n) :
    return -add(moebius(d)*log(n/d) for d in prime_divisors(n))
[round(lambda_prime(n),4) for n in (1..19)]

0.0, 0.0, 0.0, 0.6931, 0.0, 1.7918, 0.0, 1.3863, 1.0986, 2.3026,
0.0, 3.1781, 0.0, 2.6391, 2.7081, 2.0794, 0.0, 3.989, 0.0, 

Applying the exponential we get this time:

1, 1, 1, 2, 1, 6, 1, 4, 3, 10, 1, 24, 1, 14, 15, 8, 1, 54, 1, .. 

This sequence is not in the database; so we will define a new sequence.

def A205959(n) :
    return simplify(exp(
        -add(moebius(d)*log(n/d) for d in prime_divisors(n))))

Again the many 1's in this sequence are striking. The sequence of their positions is A008578 and bears the lovely name "Prime numbers at the beginning of the 20th century".  ;-))

Since A205959 is such a nice sequence we introduce a special notation for it.

Next we look at n / exp Λp(n).

[n / A205959(n) for n in (1..36)]

1, 2, 3, 2, 5, 1, 7, 2, 3, 1, 11, 1/2, 13, 1, 1, 2, 17, 1/3, 19,
1/2, 1, 1, 23, 1/4, 5, 1, 3, 1/2, 29, 1/30, 31, 2, 1, 1, 1, 1/6,

At what indices does exp Λp(n) not divide n?

12, 18, 20, 24, 28, 30, 36, 40, 42, 44, 45, 48, 50, 52, 54, 56

Looking it up we find a sequence by Reinhard Zumkeller, A102467. It is called "numbers such that the sum of numbers of prime factors with and without repetitions does not equal the number of divisors. " According to a comment of M. F. Hasler these "are the numbers which are neither prime powers (>1) nor semiprimes." (See also Hasler's remark in A102466.)

By the way, Sage is well equipped to implement this sequence using the build in 'sloane.Axxxxxx' functions. Another reason to switch to Sage.

def is_A102467(n) :
    return bool(sloane.A001221(n) <> 1 and sloane.A001222(n) <> 2)

def A102467_list(n) :
    return [k for k in (1..n) if is_A102467(k)]

Another natural question is: What are the fixed points of exp Λp(n)? In other words, for what n is exp Λp(n) = n?

1, 6, 10, 14, 15, 21, 22, 26, 33, 34, 35, 38, 39, 46, 51, 55,... 

These are the numbers that are the product of two distinct primes, A006881 (disregarding the 1).

Next let us look at the summatory function of Λp(n). It is well known that the summatory function of the von Mangoldt function is the Chebyshev function

and exp ψ(n) is lcm{1,2,3,...n}, A003418. Similarly let

What is exp ψp(n)? Perhaps the reader wishes to figure it out by own work. The answer is not A025527.

The von Mangoldt transformation

In the section 'transformations based on connection coefficients' we looked at sums which had as blueprint

    T(s, n) = Sum( CC(n, k)*s(k), k in (1..n) ).

The connection coefficients CC(n, k) were fixed and the sequence s(k) was given as input. Now let us turn the tables. Let us fix the one-parameter sequence and supply sequences with two parameters as input (we might think of  'integer triangles').

We better look for sequences with a strong arithmetical power as kernels for these transformations to make the results insightful. Luckily enough, we just have found such a function: von Mangoldt's function certainly qualifies for playing the rôle of the fixed sequence. So let's try this.

def Lambda(n) :
    return add(moebius(d)*log(n/d) for d in divisors(n))  

def mangoldt_trans(T, n) :
    return simplify(exp(add(Lambda(k)*T(n,k) for k in (1..n))))

Next, what input should we provide? In this regard we just learned a lesson while contemplating the Euler transformation: the most simple sequences, the 01-sequences delivered the most interesting results, euler_trans ([1,0,1,0,1,0,1,..]), euler_trans ([0,1,0,1,0,1,0,..]). And via Gordon's theorem the sequences [1,0,1,1,0,1,0..], [1,1,0,0,1,1,0..], etc. allowed us to catch a glimps of the celebrated Rogers-Ramanujan identities. So let us proceed on this purely formal route.

There certainly are many 01-sequences already in the database. However it seems that the compartment of 01-triangles never had been developed systematically. So let us first add some of them to the OEIS.

def even_quot(n,k)  : return int(is_even(n//k))
def odd_quot(n,k)   : return int(is_odd(n//k))
def is_divisor(k,n) : return int((n % k) == 0)
def is_primedivisor(k,n) : return int(is_prime(k) and (n%k)==0)

And indeed we find some interesting transforms.

Boolean triangle Mangoldt transform
even quot   A180000 LCM/swing
odd quot   A056040 swing
is divisor   A014963 exp mangoldt
is prime divisor   A089026 primes indicator
iota(n, k) = 1   A003418 LCM
mu(n, k) = -mu(k)   A034386 primorial

The LCM(n) and it's prime variant

First a remark on notation. We will write LCM(n) for lcm(1,2,3,...n). Here LCM is the name of a mathematical function and should not be confused with the acronym for the least common multiple. The LCM of n is n, LCM(n) is 1,2,6,12,... for n=1,2,3,4,...

The von Mangoldt transformation obviously has a prime variant.

def Lambda_prime(n) :
    return -add(moebius(d)*log(n/d) for d in prime_divisors(n))

def mangoldt_primetrans(T, n) :
    return simplify(exp(add(Lambda_prime(k)*T(n,k) for k in (1..n))))

By the iota-triangle we understand the triangle which consists completely out of 1s. The iota-triangle is mapped by the Mangoldt transformation on the LCM. Accordingly we will call the value of the iota-triangle under the Mangoldt prime transformation the prime LCM and denote it by LCMp.

def lcm_prime(n) :
    return mangoldt_primetrans(iota, n)

Making the definition explicit we can write it as

def lcm_prime(n) :
    return simplify(exp(-add(add(moebius(p)*log(k/p) 
        for p in prime_divisors(k)) for k in (1..n))))    

1, 1, 1, 1, 2, 2, 12, 12, 48, 144, 1440, 1440, 34560, 34560, 483840,
7257600, 58060800, 58060800, 3135283200, 3135283200, 125411328000, ..

This is sequence A205957. Comparing this with the prime variant of  Chebyshev's psi function as defined above we can finally state

These are also for n > 0 the partial products of A205959

1, 1, 1, 2, 1, 6, 1, 4, 3, 10, 1, 24, 1, 14, 15,...

A216152 gives the distinct values of A205957

1, 2, 12, 48, 144, 1440, 34560, 483840, 7257600,...

which in turn are the partial products of A216153

1, 2, 6, 4, 3, 10, 24, 14, 15, 8, 54, 40, 21, 22, 96,...
Let a(n) A216153(n) and b(n) = A018252(n) the nonprime numbers.
If n = 1 then 
   * a(n) = b(n) = 1;
else if a(n) < b(n) then
   * a(n) is a cototient of consecutive pure powers of primes (A053211),
   * b(n) is a prime power with exponent > 1 (A025475),
   * b(n)/a(n) is a prime root of n-th nontrivial prime power (A025476);
else if a(n) > b(n) then
   * b(n) is a number which is neither a prime power nor 
     a semiprime (A102467);
else if a(n) = b(n) then
   * a(n) is the product of two distinct primes (A006881).

A216152 are also the A205957(n) where n is a nonprime number. Moreover, the first few A216152(n) are highly totient numbers (A097942) and products of distinct factorials (A058295). Is this true for arbitrary A216152(n)?