# Dirichlet Convolution

The Dirichlet convolution is a special kind of convolution that appears as a very useful tool in number theory. It operates on the set of arithmetic functions.

### Challenge

Given two arithmetic functions $$\f,g\$$ (i.e. functions $$\f,g: \mathbb N \to \mathbb R\$$) compute the Dirichlet convolution $$\(f * g): \mathbb N \to \mathbb R\$$ as defined below.

### Details

• We use the convention $$\ 0 \notin \mathbb N = \{1,2,3,\ldots \}\$$.
• The Dirichlet convolution $$\f*g\$$ of two arithmetic functions $$\f,g\$$ is again an arithmetic function, and it is defined as $$(f * g)(n) = \sum_\limits{d|n} f\left(\frac{n}{d}\right)\cdot g(d) = \sum_{i\cdot j = n} f(i)\cdot g(j).$$ (Both sums are equivalent. The expression $$\d|n\$$ means $$\d \in \mathbb N\$$ divides $$\n\$$, therefore the summation is over the natural divisors of $$\n\$$. Similarly we can subsitute $$\ i = \frac{n}{d} \in \mathbb N, j =d \in \mathbb N \$$ and we get the second equivalent formulation. If you're not used to this notation there is a step by step example at below.) Just to elaborate (this is not directly relevant for this challenge): The definition comes from computing the product of Dirichlet series: $$\left(\sum_{n\in\mathbb N}\frac{f(n)}{n^s}\right)\cdot \left(\sum_{n\in\mathbb N}\frac{g(n)}{n^s}\right) = \sum_{n\in\mathbb N}\frac{(f * g)(n)}{n^s}$$
• The input is given as two black box functions. Alternatively, you could also use an infinite list, a generator, a stream or something similar that could produce an unlimited number of values.
• There are two output methods: Either a function $$\f*g\$$ is returned, or alternatively you can take take an additional input $$\n \in \mathbb N\$$ and return $$\(f*g)(n)\$$ directly.
• For simplicity you can assume that every element of $$\ \mathbb N\$$ can be represented with e.g. a positive 32-bit int.
• For simplicity you can also assume that every entry $$\ \mathbb R \$$ can be represented by e.g. a single real floating point number.

### Examples

Let us first define a few functions. Note that the list of numbers below each definition represents the first few values of that function.

• the multiplicative identity (A000007) $$\epsilon(n) = \begin{cases}1 & n=1 \\ 0 & n>1 \end{cases}$$ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
• the constant unit function (A000012)$$\mathbb 1(n) = 1 \: \forall n$$ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
• the identity function (A000027) $$id(n) = n \: \forall n$$ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, ...
• the Möbius function (A008683) $$\mu(n) = \begin{cases} (-1)^k & \text{ if } n \text{ is squarefree and } k \text{ is the number of Primefactors of } n \\ 0 & \text{ otherwise } \end{cases}$$ 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, ...
• the Euler totient function (A000010) $$\varphi(n) = n\prod_{p|n} \left( 1 - \frac{1}{p}\right)$$ 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8, ...
• the Liouville function (A008836) $$\lambda (n) = (-1)^k$$ where $$\k\$$ is the number of prime factors of $$\n\$$ counted with multiplicity 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, ...
• the divisor sum function (A000203) $$\sigma(n) = \sum_{d | n} d$$ 1, 3, 4, 7, 6, 12, 8, 15, 13, 18, 12, 28, 14, 24, 24, 31, 18, 39, 20, ...
• the divisor counting function (A000005) $$\tau(n) = \sum_{d | n} 1$$ 1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6, 2, 4, 4, 5, 2, 6, 2, 6, 4, 4, 2, 8, ...
• the characteristic function of square numbers (A010052) $$sq(n) = \begin{cases} 1 & \text{ if } n \text{ is a square number} \\ 0 & \text{otherwise}\end{cases}$$ 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, ...

Then we have following examples:

• $$\ \epsilon = \mathbb 1 * \mu \$$
• $$\ f = \epsilon * f \: \forall f \$$
• $$\ \epsilon = \lambda * \vert \mu \vert \$$
• $$\ \sigma = \varphi * \tau \$$
• $$\ id = \sigma * \mu\$$ and $$\ \sigma = id * \mathbb 1\$$
• $$\ sq = \lambda * \mathbb 1 \$$ and $$\ \lambda = \mu * sq\$$
• $$\ \tau = \mathbb 1 * \mathbb 1\$$ and $$\ \mathbb 1 = \tau * \mu \$$
• $$\ id = \varphi * \mathbb 1 \$$ and $$\ \varphi = id * \mu \$$

The last for are a consequence of the Möbius inversion: For any $$\f,g\$$ the equation $$\ g = f * 1\$$ is equivalent to $$\f = g * \mu \$$.

### Step by Step Example

This is an example that is computed step by step for those not familiar with the notation used in the definition. Consider the functions $$\f = \mu\$$ and $$\g = \sigma\$$. We will now evaluate their convolution $$\\mu * \sigma\$$ at $$\ n=12\$$. Their first few terms are listed in the table below.

$$\begin{array}{c|ccccccccccccc} f & f(1) & f(2) & f(3) & f(4) & f(5) & f(6) & f(7) & f(8) & f(9) & f(10) & f(11) & f(12) \\ \hline \mu & 1 & -1 & -1 & 0 & -1 & 1 & -1 & 0 & 0 & 1 & -1 & 0 \\ \sigma & 1 & 3 & 4 & 7 & 6 & 12 & 8 & 15 & 13 & 18 & 12 & 28 \\ \end{array}$$

The sum iterates over all natural numbers $$\ d \in \mathbb N\$$ that divide $$\n=12\$$, thus $$\d\$$ assumes all the natural divisors of $$\n=12 = 2^2\cdot 3\$$. These are $$\d =1,2,3,4,6,12\$$. In each summand, we evaluate $$\g= \sigma\$$ at $$\d\$$ and multiply it with $$\f = \mu\$$ evaluated at $$\\frac{n}{d}\$$. Now we can conclude

$$\begin{array}{rlccccc} (\mu * \sigma)(12) &= \mu(12)\sigma(1) &+\mu(6)\sigma(2) &+\mu(4)\sigma(3) &+\mu(3)\sigma(4) &+\mu(2)\sigma(6) &+\mu(1)\sigma(12) \\ &= 0\cdot 1 &+ 1\cdot 3 &+ 0 \cdot 4 &+ (-1)\cdot 7 &+ (-1) \cdot 12 &+ 1 \cdot 28 \\ &= 0 & + 3 & 1 0 & -7 & - 12 & + 28 \\ &= 12 \\ & = id(12) \end{array}$$

# Lean, 1081009578 75 bytes

def d(f g:_->int)(n):=(list.iota n).foldr(λd s,ite(n%d=0)(s+f d*g(n/d))s)0


Try it online!

More testcases with all of the functions.

• is lambda really more expensive than four bytes for fun ? Nov 16, 2018 at 21:50
• lambda is three bytes, I suppose Nov 16, 2018 at 21:50
• I think it's two in UTF8 (greek is pretty low unicode) Nov 16, 2018 at 21:53
• You're right. I also golfed the import Nov 16, 2018 at 21:57
• I also used cond to save 5 bytes Nov 16, 2018 at 22:00

(f!g)n=sum[f i*g(div n i)|i<-[1..n],mod n i<1]


Try it online!

Thanks to flawr for -6 bytes and a great challenge! And thanks to H.PWiz for another -6!

• Simpler is shorter here Nov 17, 2018 at 3:04
• @H.PWiz That's pretty clever - I didn't even think of doing it that way!
– user45941
Nov 17, 2018 at 4:19

# Wolfram Language (Mathematica), 17 bytes

Of course Mathematica has a built-in. It also happens to know many of the example functions. I've included some working examples.

DirichletConvolve


Try it online!

# Python 3, 59 bytes

lambda f,g,n:sum(f(d)*g(n//d)for d in range(1,n+1)if 1>n%d)


Try it online!

• Is // really needed instead of /? Nov 16, 2018 at 20:08
• / would produce floats right? Nov 16, 2018 at 20:09
• Because d is a divisor of n by definition, the fractional part of n/d is zero, so there shouldn't be any issues with floating point arithmetic. Floats with fractional part zero are close enough to ints for Pythonic purposes, and the output of the function is a real number, so doing n/d instead of n//d should be fine.
– user45941
Nov 17, 2018 at 4:28

D,g,@~,$z€¦~¦* D,f,@@@,@b[VdF#B]dbRzGb]$dbL$@*z€g¦+  Try it online! Takes two pre-defined functions as arguments, plus $$\n\$$, and outputs $$\(f * g)(n)\$$ ## How it works D,g, ; Define a helper function,$g
@~,	; $g takes a single argument, an array, and splats that array to the stack ;$g takes the argument e.g. [[τ(x) φ(x)] [3 4]]
; STACK : 			[[τ(x) φ(x)] [3 4]]
$z ; Swap and zip: [[3 τ(x)] [4 φ(x)]] €¦~ ; Reduce each by execution: [[τ(3) φ(4)]] ¦* ; Take the product and return: τ(3)⋅φ(4) = 4 D,f, ; Define the main function,$f
@@@,	; $f takes three arguments: φ(x), τ(x) and n (Let n = 12) ; STACK: [φ(x) τ(x) 12] @ ; Reverse the stack: [12 τ(x) φ(x)] b[V ; Pair and save: [12] Saved: [τ(x) φ(x)] dF#B] ; List of factors: [[1 2 3 4 6 12]] dbR ; Copy and reverse: [[1 2 3 4 6 12] [12 6 4 3 2 1]] z ; Zip together: [[[1 12] [2 6] [3 4] [4 3] [6 2] [12 1]]] Gb] ; Push Saved: [[[1 12] [2 6] [3 4] [4 3] [6 2] [12 1]] [[τ(x) φ(x)]]]$dbL	; Number of dividors:		[[[τ(x) φ(x)]] [[1 12] [2 6] [3 4] [4 3] [6 2] [12 1]] 6]
$@* ; Repeat: [[[1 12] [2 6] [3 4] [4 3] [6 2] [12 1]] [[τ(x) φ(x)] [τ(x) φ(x)] [τ(x) φ(x)] [τ(x) φ(x)] [τ(x) φ(x)] [τ(x) φ(x)]]] z ; Zip: [[[τ(x) φ(x)] [1 12]] [[τ(x) φ(x)] [2 6]] [[τ(x) φ(x)] [3 4]] [[τ(x) φ(x)] [4 3]] [[τ(x) φ(x)] [6 2]] [[τ(x) φ(x)] [12 1]]] €g ; Run$g over each subarray:	[[4 4 4 6 4 6]]
¦+	; Take the sum and return:	28


# R, 58 bytes

function(n,f,g){for(i in (1:n)[!n%%1:n])F=F+f(i)*g(n/i)
F}


Try it online!

Takes n, f, and g. Luckily the numbers package has quite a few of the functions implemented already.

If vectorized versions were available, which is possible by wrapping each with Vectorize, then the following 45 byte version is possible:

# R, 45 bytes

function(n,f,g,x=1:n,i=x[!n%%x])f(i)%*%g(n/i)


Try it online!

# APL (Dyalog Classic), 20 bytes

{(⍺⍺¨∘⌽+.×⍵⍵¨)∪⍵∨⍳⍵}


with ⎕IO←1

Try it online!

Easy to solve, hard to test - generally not my type of challenge. Yet, I enjoyed this one very much!

{ } defines a dyadic operator whose operands ⍺⍺ and ⍵⍵ are the two functions being convolved; ⍵ is the numeric argument

∪⍵∨⍳⍵ are the divisors of ⍵ in ascending order, i.e. unique (∪) of the LCMs (∨) of ⍵ with all natural numbers up to it (⍳)

⍵⍵¨ apply the right operand to each

⍺⍺¨∘⌽ apply the left operand to each in reverse

+.× inner product - multiply corresponding elements and sum

The same in ngn/apl looks better because of Unicode identifiers, but takes 2 additional bytes because of 1-indexing.

• Pretty sure it takes 27 additional bytes in ngn/apl... Nov 17, 2018 at 11:42

# Jelly, 9 bytes

ÆDṚÇ€ḋÑ€Ʋ


Try it online!

Line at the top is the main line of $$\f\$$, line at the bottom is the main line of $$\g\$$. $$\n\$$ is passed as an argument to this function.

# Swift 4,  74 70  54 bytes

{n in(1...n).map{n%$0<1 ?f(n/$0)*g(\$0):0}.reduce(0,+)}


Try it online!

# JavaScript (ES6), 47 bytes

Takes input as (f)(g)(n).

f=>g=>h=(n,d=n)=>d&&!(n%d)*f(n/d)*g(d)+h(n,d-1)


Try it online!

### Examples

liouville =
n => (-1) ** (D = (n, k = 2) => k > n ? 0 : (n % k ? D(n, k + 1) : 1 + D(n / k, k)))(n)

mobius =
n => (M = (n, k = 1) => n % ++k ? k > n || M(n, k) : n / k % k && -M(n / k, k))(n)

sq =
n => +!((n ** 0.5) % 1)

identity =
n => 1

// sq = liouville * identity
console.log([...Array(25)].map((_, n) => F(liouville)(identity)(n + 1)))

// liouville = mobius * sq
console.log([...Array(20)].map((_, n) => F(mobius)(sq)(n + 1)))


# C (gcc), 108 bytes

#define F float
F c(F(*f)(int),F(*g)(int),int n){F s=0;for(int d=0;d++<n;)if(n%d<1)s+=f(n/d)*g(d);return s;}


Straightforward implementation, shamelessly stolen from Leaky Nun's Python answer.

Ungolfed:

float c(float (*f)(int), float (*g)(int), int n) {
float s = 0;
for(int d = 1; d <= n;++d) {
if(n % d == 0) {
s += f(n / d) * g(d);
}
}
return s;
}


Try it online!

• 95 Jan 5, 2020 at 17:12

## F#, 72 bytes

let x f g n=Seq.filter(fun d->n%d=0){1..n}|>Seq.sumBy(fun d->f(n/d)*g d)


Takes the two functions f and g and a natural number n. Filters out the values of d that do not naturally divide into n. Then evaluates f(n/d) and g(d), multiples them together, and sums the results.

# Pari/GP, 32 bytes

(f,g,n)->sumdiv(n,d,f(n/d)*g(d))


Try it online!

There is a built-in dirmul function, but it only supports finite sequences.

# APL(NARS), 47 chars, 94 bytes

{(m⍵[1])×n⍵[2]}{+/⍺⍺¨{k←⍳⍵⋄(⍵÷b),¨b←k/⍨0=k∣⍵}⍵}


where m and n are the function one has to use (this because i don't know how to call one array of function in one function in APL). Using the example above the multiplication of Mobius function (here it is 12π) and sum of divisors function (here it is 11π) for value 12 that multiplication would be:

  {(12π⍵[1])×11π⍵[2]}{+/⍺⍺¨{k←⍳⍵⋄(⍵÷b),¨b←k/⍨0=k∣⍵}⍵}12
12


if one has to calculate some other value:

  {(12π⍵[1])×11π⍵[2]}{+/⍺⍺¨{k←⍳⍵⋄(⍵÷b),¨b←k/⍨0=k∣⍵}⍵}1002
1002
{(12π⍵[1])×11π⍵[2]}{+/⍺⍺¨{k←⍳⍵⋄(⍵÷b),¨b←k/⍨0=k∣⍵}⍵}1001
1001
{(12π⍵[1])×11π⍵[2]}{+/⍺⍺¨{k←⍳⍵⋄(⍵÷b),¨b←k/⍨0=k∣⍵}⍵}20000x
20000


One can see if for example the first 2000 number the function result is the identity

  (⍳2000)≡{(12π⍵[1])×11π⍵[2]}{+/⍺⍺¨{k←⍳⍵⋄(⍵÷b),¨b←k/⍨0=k∣⍵}⍵}¨⍳2000
1