7
\$\begingroup\$

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]
\$\endgroup\$
  • \$\begingroup\$ I think solving a linear system with inbuilt functions should be prohibited completely. \$\endgroup\$ – Optimizer Nov 29 '14 at 6:39
  • \$\begingroup\$ @Optimizer Using built-ins to solve linear systems is already banned; what more do you mean? Prohibiting inverting matrices? In any case, I'm not going to make changes at this point. I agree that matrix-based solutions are hard to beat, which I hadn't foreseen, so please consider it a separate league if it makes you feel better. \$\endgroup\$ – xnor Nov 29 '14 at 8:59
  • \$\begingroup\$ See the J answer. Aolving Ax=B for a matrix is equivalent to solving system of equations only. \$\endgroup\$ – Optimizer Nov 29 '14 at 20:52
5
\$\begingroup\$

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
\$\endgroup\$
  • \$\begingroup\$ Sorry, but you're not allowed to use methods that solve linear systems. You can perform matrix operations though (or is that what your code is already doing?). \$\endgroup\$ – xnor Nov 29 '14 at 8:54
  • \$\begingroup\$ @xnor Yes, it's only matrix inversion. It states in the the documentation that b%.A is shorthand for (%.A) mp b where mp is the matrix product, and %.A is matrix inversion on A. \$\endgroup\$ – algorithmshark Nov 29 '14 at 14:18
2
\$\begingroup\$

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.

\$\endgroup\$
2
\$\begingroup\$

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]
\$\endgroup\$
2
\$\begingroup\$

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]
\$\endgroup\$
2
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ How is this solving the system of equations? If it's just using the built-in CAS solver, that's not allowed, but if you're inverting the matrix, that's fine. My apologies that I edited in that restriction after posting the question, so you might not have seen it. \$\endgroup\$ – xnor Nov 28 '14 at 23:48
  • \$\begingroup\$ It is using a built in numerical solver (you can assume it is the gauss algorithm), Matlab/Octave are numerical computing languages that do not allow algebraic manipulations - so is this still ok? \$\endgroup\$ – flawr Nov 28 '14 at 23:49
  • \$\begingroup\$ I think I need to understand what (i/~mod(x,x') is doing. Could you please explain? \$\endgroup\$ – xnor Nov 28 '14 at 23:52
  • \$\begingroup\$ ´i´ is the input vector (as a row vector). / is divides the vector by the matrix and ~mod(x,x') represents the matrix. (I'll write more in my answer, just noticed the comment is too short.) \$\endgroup\$ – flawr Nov 28 '14 at 23:54
  • \$\begingroup\$ OK, that's totally fine. The code expresses an explicit solution to the matrix equation; it's not being implicitly solved. \$\endgroup\$ – xnor Nov 28 '14 at 23:56
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.