How can we compute the multiplicative partition function

Solution 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;

Solution 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.$