Mobius inversion

Question

Your goal is to implement Möbius inversion, a linear operation on a sequence of numbers that is used in combinatorics and number theory. Fewest bytes wins.

What is Möbius inversion?

Given a (finite) integer sequence f_1, f_2, ..., f_N, we can compute its divisor-partial-sum g, whose nth element g_n is the sum of all terms f_d whose index d is a divisor of n.

g_n = sum_(d with d|n) [f_d]

g_1 = f_1
g_2 = f_1 + f_2
g_3 = f_1 + f_3
g_4 = f_1 + f_2 + f_ 4
g_5 = f_1 + f_5
g_6 = f_1 + f_2 + f_3 + f_6
...

The sequence g has the same length as f.

Möbius inversion is the process of inverting this operation -- recovering f from g. This inverse is unique and produces an integer sequence. A explicit but complicated way to compute it is the Möbius inversion formula, which expresses f_n as a weighted sum over divisors d of n:

f_n = sum_(d with d|n) [μ(d) * g_(n/d)]

and μ is the Möbius function:

μ(n) = 1 if n is a square-free positive integer with an even number of prime factors.
μ(n) = −1 if n is a square-free positive integer with an odd number of prime factors.
μ(n) = 0 if n has a squared prime factor.

Program requirements

Write a function or program to perform Möbius inversion. Fewest bytes wins. Built-in methods to perform Möbius inversion, do sequence convolution, compute the Möbius function (or related thing like checking squarefreeness), or solve systems of equations are not allowed.

Input

A sequence of integers of length at least 1. Can be any format that represents an ordered sequence: array, list, stack, token-separated string, etc. Unordered maps from the index to the element like dictionaries and functions are not OK. You're not given the length of the input sequence separately.

Output

Return or print the output sequence in any format, with same restrictions as input. It should be the same length as the input sequence.

Test cases

Input: [1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6, 2, 4, 4, 5, 2, 6, 2, 6]
Output:[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

Input: [0, -3, 6, 5, 1, 7, 5, 0, 4, 0]
Output: [0, -3, 6, 8, 1, 4, 5, -5, -2, 2]

Input: [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Output: [1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1]

Answer 1

J - 13 char

Anonymous function Mobius-inverting a sequence.

(%.0=#\|~/#\)

Here's what's going on:

• #\ gets indices from 1 to the length of the input list, inclusive.
• |~/ makes a table of indices modulo themselves.

• 0= checks for equality to zero, 1 for true, 0 for false. On the table of moduli, this checks divisibility, returning a matrix like

  1 0 0 0 0 0 0 0
  1 1 0 0 0 0 0 0
  1 0 1 0 0 0 0 0
  1 1 0 1 0 0 0 0
  1 0 0 0 1 0 0 0
  1 1 1 0 0 1 0 0
  1 0 0 0 0 0 1 0
  1 1 0 1 0 0 0 1
  

• Finally, %. takes the input sequence and this matrix, and solves it as a linear system: x equal to b %. A is the solution to Ax=b if it exists. Our b is the input sequence and A is this divisibility matrix above.

In use: (Note that J spells its negative sign _)

   (%.0=#\|~/#\) 1 2 2 3 2 4 2 4 3 4 2 6 2 4 4 5 2 6 2 6
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
   (%.0=#\|~/#\) 0 _3 6 5 1 7 5 0 4 0
0 _3 6 8 1 4 5 _5 _2 2
   (%.0=#\|~/#\) 1 0 0 0 0 0 0 0 0 0 0 0 0
1 _1 _1 0 _1 1 _1 0 0 1 _1 0 _1

Answer 2

Python, 99 chars

def I(g):
 f,e=[],enumerate
 for n,v in e(g):f+=[v-sum(x for i,x in e(f)if(n+1)%(i+1)<1)]
 return f

Computes f from g incrementally, so we can use the fs we've already calculated to compute larger fs.

Answer 3

CJam, 67 bytes

This can be golfed further, but given that so many languages can solve linear equations with 1 character, I don't see a point in answering with any non-maths oriented language.

l~:F,,{):N,:){N\%!},{__mF{W=1>},,{;0}{_mF,\1=m2%W1?}?N@/(F=*}%:+}%p

This is the exact implementation of the Möbius inversion formula.

Try it online here

Takes input from STDIN in the following format:

[1 0 0 0 0 0 0 0 0 0 0 0 0]

and gives output line

[1 -1 -1 0 -1 1 -1 0 0 1 -1 0 -1]

Answer 4

Haskell, 61

h g=[v-sum[h g!!n|n<-[0..i-2],rem i(n+1)<1]|(i,v)<-zip[1..]g]

this solves the equations one at a time using the already known values of f, using laziness.

this is inefficient in space and time complexity because it recomputes the result several times ( it's like the difference between fix f=f(fix f) and fix f=r where r=f r). a simple fix:

h g=r where [v-sum[r!!n|n<-[0..i-2],rem i(n+1)<1]|(i,v)<-zip[1..]g]

Answer 5

Matlab/Octave (56)(53)(51)

51 chars:

i=input('');x=1:numel(i);disp(i/~bsxfun(@mod,x,x'))

53 chars:

i=input('');x=meshgrid(1:numel(i));disp(i/~mod(x,x'))

I basically set up system of linear equations in matrix form where the matrix looks like this (depending on the length, here the length is 6:

1 0 0 0 0 0 
1 1 0 0 0 0 
1 0 1 0 0 0 
1 1 0 1 0 0 
1 0 0 0 1 0 
1 1 1 0 0 1

Then I just solve the corresponding system of equations.

The command meshgrid(x) is primarly inteded for producing matrices like

1 1 1 1
2 2 2 2
3 3 3 3
4 4 4 4

and

1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4

which can then be used as e.g. x and y coordinates for plotting 2d functions. I just use one of them for ceating a matrix with where the k-th column has is 0 every k-th element via elementwise remainder division (modulo). ~ as the logical not which replaces nonzero values with zeros and vice versa, which then results in the matrix above.

Answer 6

Python 2, 84 chars

g=input()
i=1
for _ in g:g[i-1]-=sum(g[j-1]*(i%j<1)for j in range(1,i));i+=1
print g

Iterates through the entries, making each one what it needs to be for the sum to work out. Modifies the input list directly rather than making a new list.