How can we compute the multiplicative partition function

Question

Please how can we compute the multiplicative partition function. For example $24$, has precisely $6$ valid factorizations: $2\cdot2\cdot2\cdot3$, $2\cdot2\cdot6$, $2\cdot3\cdot4$, $2\cdot12$, $3\cdot8$, and $4\cdot6$.

Answer 1

There is a formula that depends only on the shape of the prime factorization of $n$, i.e. the exponents of the primes in the factorization. This formula is obtained with the Polya Enumeration Theorem (PET) applied to a certain repertoire. It includes the factorization into a product of one term, namely $n$.

Let $$n= p_1^{v_1} p_2^{v_2} \times \cdots \times p_q^{v_q}$$ be the prime factorization of $n$ and recall the function $\Omega(n)$ which is given by $$\Omega(n) = \sum_{k=1}^q v_k,$$ i.e. it counts the prime factors according to their multiplicities.

The repertoire that goes into PET is given by $$Q = -1 + \prod_{k=1}^q \sum_{v=0}^{v_k} u_{p_k}^v = -1 + \prod_{k=1}^q \frac{u_{p_k}^{v_k+1}-1}{u_{p_k}-1}$$ which is an ordinary generating function in the set of variables $\{u_{p_k}\}.$

With these settings the value $q(n)$ of the multiplicative partition function is given by $$q(n) = [u_{p_1}^{v_1} u_{p_2}^{v_2}\cdots u_{p_q}^{v_q}] \sum_{m=1}^{\Omega(n)} Z(S_m)(Q),$$ where $Z(S_m)$ is the cycle index of the symmetric group acting on $m$ elements, which is computed from $$Z(S_0) = 1 \quad \text{and}\quad Z(S_m) = \frac{1}{m} \sum_{l=1}^m a_l Z(S_{m-l}).$$

This produces the following sequence starting at $n=2$ and extending to $n=64$ $$1, 1, 2, 1, 2, 1, 3, 2, 2, 1, 4, 1, 2, 2, 5, 1, 4, 1, 4,\\ 2, 2, 1, 7, 2, 2, 3, 4, 1, 5, 1, 7, 2, 2, 2, 9, 1, 2, 2, 7,\\ 1, 5, 1, 4, 4, 2, 1, 12, 2, 4, 2, 4, 1, 7, 2, 7, 2, 2, 1, 11, 1, 2, 4, 11,\ldots$$ which points to OEIS A001055.

The complexity of the repertoire that goes into the cycle index is $\tau(n).$ Remark as of Tue Jan 7 00:40:27 CET 2014. This algorithm is definitely not as good as what was posted at the Math Overflow link. Please consult my second answer for a fast algorithm.

Here is the Maple code that was used to compute the above values:

with(numtheory);
with(group):
with(combinat):

pet_cycleind_symm :=
proc(n)
        option remember;

        if n=0 then return 1; fi;

        expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

pet_varinto_cind :=
proc(poly, ind)
           local subs1, subs2, polyvars, indvars, v, pot, res;

           res := ind;

           polyvars := indets(poly);
           indvars := indets(ind);

           for v in indvars do
               pot := op(1, v);

               subs1 :=
               [seq(polyvars[k]=polyvars[k]^pot,
               k=1..nops(polyvars))];

               subs2 := [v=subs(subs1, poly)];

               res := subs(subs2, res);
           od;

           res;
end;

v :=
proc(n)
        option remember;
        local fact, rep, gf, fcount, p, res;

        fact := op(2, ifactors(n));

        rep := mul(add(u[fact[k][1]]^p, p=0..fact[k][2]),
                   k=1..nops(fact));
        rep := rep-1;

        res := 0;
        for fcount to bigomega(n) do
            gf := pet_varinto_cind(rep, pet_cycleind_symm(fcount));
            gf := expand(gf);

            for p to nops(fact) do
                gf := coeff(gf, u[fact[p][1]], fact[p][2]);
            od;

            res := res+gf;
        od;

        res;
end;

Answer 2

For the purpose of completing this discussion and providing an effective solution I am sending a dynamic-programming implementation in Maple and in Perl. This is the Maple code. Addendum Tue Jan 7 08:05:14 CET 2014. Further optimization is possible as this statistic only depends on the shape of the factorization and not the primes, so the shape can be our key for use in memoization. The reader is invited to try this. Of course in this setup you only profit from memoization if you invoke the function at the main interface multiple times, so this is not a DP algorithm in the classic sense.

mp :=
proc(n)
    option remember;
    local rec, onef, res, f, asg, targ, t;

    if n=1 then return {}; fi;

    onef := op(2, ifactors(n))[1];

    rec := mp(n/onef[1]);
    if nops(rec) = 0 then return {[n=1]} fi;

    res := {};
    for f in rec do
        t := table(f);

        targ := onef[1];
        if type(t[targ], integer) then
            t[targ] := t[targ] + 1
        else
            t[targ] := 1;
        fi;

        res := res union {sort(op(op(t)))};

        for asg in f do
            t := table(f);

            targ := onef[1]*op(1, asg);
            if type(t[targ], integer) then
                t[targ] := t[targ]+1;
            else
                t[targ] := 1;
            fi;

            if t[op(1,asg)]>1 then
                t[op(1,asg)] := t[op(1,asg)]-1;
            else
                t[op(1,asg)] := evaln(t[op(1,asg)]); 
            fi;

            res := res union {sort(op(op(t)))};
        od;
    od;

    res;
end;

It can be used to compute the number of factorizations approximately up to about $9!$, after which it slows dramatically for some reason.

To get around this problem I am also sending a Perl program that implements the DP algorithm. This has the drawback that the factorization method is naive and slow to detect large primes but it illustrates the concept especially for values with many factors that have numerous factorizations e.g. for $10^{10}$ which gives $59521$ factorizations.

#! /usr/bin/perl -w
#

sub onefact {
    my ($n) = @_;

    return 2 if $n % 2 == 0;

    my $f = 3;

    while($n % $f > 0){
        $f += 2;
    }

    return $f;
}

my %memo = ( 1 => {} );

sub mp {
    my ($n) = @_;
    return $memo{$n} if exists $memo{$n};

    my $res = {};

    my $onef = onefact($n);
    if($n/$onef == 1){
        $res->{"$n=1"} = 1;
    }
    else{
        sub tbl2str {
            my ($tref) = @_;
            my @keyl = map {
                $_ . "=" . $tref->{$_}
            } sort(keys(%$tref));

            return join('-', @keyl);
        }

        foreach my $f (keys %{mp($n/$onef)}){
            my @asigs = map {
                split(/=/);
            } split(/-/, $f);

            my (%tblsrc) = @asigs;

            my (%tbl) = @asigs;
            if(exists($tbl{$onef})){
                $tbl{$onef} = $tbl{$onef} + 1;
            }
            else{
                $tbl{$onef} = 1;
            }

            $res->{tbl2str(\%tbl)} = 1;

            for my $a (keys %tblsrc){
                %tbl = %tblsrc;

                my $targ = $onef*$a;
                if(exists($tbl{$targ})){
                    $tbl{$targ}++;
                }
                else{
                    $tbl{$targ} = 1;
                }

                if($tbl{$a}>1){
                    $tbl{$a}--;
                }
                else{
                    delete $tbl{$a};
                }

                $res->{tbl2str(\%tbl)} = 1;
            }
        }
    }

    $memo{$n} = $res;
    return $res;
}


MAIN: {
    while(my $n = shift){
        my $res = scalar(keys %{mp($n)});
        my $msize = scalar(keys %memo);
        print "$n $res [$msize]\n";
    }
}

This program is very fast and can be used to compute the number of factorizations even of large factorials, e.g. for $8!, 9!, 10!$ and $11!$ it computes in a matter of seconds the values $$2116,11830, 70520, 425240,2787810\ldots$$ which points us to OEIS A076716. Perl took about eight minutes to compute the value for $12!.$ In spite of the naive factorization the Perl program can be used on large primes or products of them, e.g. try $104729$ or $62773913 = 7919\times 7927.$

Answer 3

I feel compelled to add another answer as the code from my initial effort is really quite poor and I am not only referring to the grade school factorization routine. I asked the kind folks at MaplePrimes for advice and it turns out that the multiplicative partition function was the subject of an extensive discussion there and a collaborative effort led to a simple and effective algorithm being developed. The key observation in this algorithm is that we can generate all partitions by starting with the single factor partition and calculating the two-factor, three-factor etc. partitions by repeatedly splitting the largest element into two factors and admitting the resulting partition into the new collection if both factors are at least as large as the factor at the end of the partition whose last element was split. That is all and this simple observation makes it possible to compute very large sets of multiplicative partitions quickly and effectively. The credit goes to the MaplePrimes Thread of course and to the people who participated in it. What is nice about this algorithm is that it maintains the partitions and their elements in sorted order automagically.

I am sending a Perl program that implements this simple algorithm (there is a Maple implementation including extensive commentary at the link I gave). It takes three types of arguments: a factorial produces the count of the multiplicative partition values at a range of factorials up to some maximum, a range with a central dash will process an entire range of values and a single numeric argument outputs the partitions for that argument.

Here are some examples.

$ ./mp-final.pl 40-60
7, 1, 5, 1, 4, 4, 2, 1, 12, 2, 4, 2, 4, 1, 7, 2, 7, 2, 2, 1, 11
$ ./mp-final.pl 10!
0, 1, 2, 7, 21, 98, 392, 2116, 11830, 70520
$ ./mp-final.pl 100
> 100
> 5 20
> 2 50
> 10 10
> 4 25
> 2 5 10
> 2 2 25
> 4 5 5
> 2 2 5 5
9

The chief drawback of this Perl skript is that the factorization is naive and there is no builtin for computing the number theoretic function $\Omega(n)$ giving the prime count including multiplicities.

This is the code, and improvement over the first version I hope although that is not saying much.

#! /usr/bin/perl -w
#

sub onefact {
    my ($n) = @_;

    return 2 if $n % 2 == 0;

    my $f = 3;
    my $z = int sqrt $n;

    while($n % $f > 0 && $f <= $z){
        $f += 2;
    }
    return $n if $f>$z;

    return $f;
}


my %div_memo;
sub div {
    my ($n) = @_;
    return $div_memo{$n} if exists $div_memo{$n};

    my $of = onefact($n);
    if($of == $n){
        $div_memo{$n} = [1, $n];
    }
    else{
        my $divs = [];
        push @$divs, @{div($n/$of)};

        my %seen = map { $_, 1 } @$divs;

        my $mx = scalar(@$divs);
        for(my $pos=0; $pos<$mx; $pos++){
            my $targ = $divs->[$pos]*$of;
            push @$divs, $targ if not exists $seen{$targ};
        }

        $div_memo{$n} = $divs;
    }

    return $div_memo{$n};
}

my %memo = ("1" => []);

sub mp {
    my ($n) = @_;
    return $memo{"$n"} if exists $memo{"$n"};

    my @result = ([$n]); my $m = 2; 
    my ($min, $max) = (0, 0);

    while(1){
        for(my $ind = $min; $ind<=$max; $ind++){
            my $mpart = $result[$ind];

            my $spl = $#$mpart;
            my $spval = $mpart->[$spl];
            my @divs = @{ div($spval) };

            foreach my $d (@divs){
                my $q = $spval/$d;

                if($d > 1 && $d <= $q && 
                   ($spl == 0 || $d >= $mpart->[-2])){
                    my $npart = [@$mpart];

                    pop @$npart;
                    push @$npart, $d, $q;

                    push @result, $npart;
                }
            }
        }

        my $lastfact = $result[-1]; $cont = undef;
        foreach my $term (@$lastfact){
            if(scalar(@{ div($term) }) != 2){
                $cont = 1;
                last;
            }
        }

        last if not defined($cont);
        $min = $max+1; $max = $#result;
        $m++;
    }

    $memo{"$n"} = \@result;
    return \@result;
}


MAIN: {
    while(my $nstr = shift){
        if(my @int = ($nstr =~ /^(\d+)-(\d+)$/)){
            for(my $n = $int[0]; $n <= $int[1]; $n++){
                print ", " if $n>$int[0];
                print scalar(@{mp($n)});
            }
        }
        elsif (my @mx = ($nstr =~ /^(\d+)!$/)){
            for(my ($n, $p) = (1, 1); $n <= $mx[0]; $n++){
                $p *= $n;

                print ", " if $p>1;
                print scalar(@{mp($p)});
            }
        }
        elsif (my @nvals = ($nstr =~ /^(\d+)$/)){
            foreach my $ent (@{mp($nvals[0])}){
                print "> @$ent\n";
            }

            print scalar(@{mp($nvals[0])});
        }
        else{
            warn "skipping argument $nstr";
        }
    }

    print "\n";
}

Answer 4

Following up on some of Marko's replies, Perl does have at least two ways to compute Ω(n), assuming you are willing to use modules: scalar Math::Prime::Util::factor, and Math::Pari::bigomega. Here is an example using modules to make simpler code. I've also added the fcnt function from A001055 that just returns the count, which lets one calculate the factorial examples faster and with reasonable amounts of memory.

use warnings;
use strict;
use Math::Prime::Util qw/divisors fordivisors is_prime/;
use Memoize;  memoize('mp');  memoize('fcnt');
$|++;

sub fcnt {
  my($n, $m) = @_;
  return 1 if $n == 1;
  my $s = 0;
  fordivisors {
    $s += fcnt($n/$_, $_) if $_ > 1 && $_ <= $m;
  } $n;
  return $s;
}

sub mp {
  my ($n) = @_;
  return [] if $n == 1;

  my @result = ([$n]);
  my ($min, $max) = (0, 0);

  while ( scalar grep { !is_prime($_) } @{$result[-1]} ) {

    for (my $ind = $min; $ind <= $max; $ind++) {
      my @mpart = @{$result[$ind]};

      my $spval = $mpart[-1];
      my $dmin = (scalar @mpart < 2) ? 2 : $mpart[-2];
      foreach my $d (divisors($spval)) {
        next unless $d >= $dmin;
        my $q = $spval/$d;
        last if $d > $q;  # Divisors are sorted
        push @result, [ @mpart[0 .. $#mpart-1], $d, $q ];
      }
    }
    ($min, $max) = ($max+1, $#result);
  }
 return \@result;
}

The MAIN section is identical other than using "print fcnt($p,$p);" in the factorial case.