Multiplicative partitions

[Charles Greathouse]

A001055(n) is the number of multiplicative partitions of n. It's marked as an easy sequence, but I can't find any program on it that can easily compute

A001055(13585283712000)

which I thought was a pretty tame example (exponents 10, 4, 3, 2, 2, 1, 1).

Is this an easy sequence that needs a better program or should it lose kw:easy? (Or maybe this is a bad example, but I don't really think it's asking that much.)

[Gareth McCaughan]

--
You received this message because you are subscribed to the Google Groups "SeqFan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to seqfan+un...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/seqfan/CANXmBjwh2508d-7rVo_B1k64_T6bmy2%3DVt-OCfwbnkyOZRf7-g%40mail.gmail.com.

My machine (desktop, 5y old but fairly capable) takes about 50s to compute A001055(13585283712000) = 4258036618 using some code rather brainlessly derived from the paper by Shamik Ghosh cited on the A001055 page. Is that "easily" enough? (I am on the fence as to whether or not it is.) I have no particular reason to think that this is anything like the optimal way to do it, but it seemed plausibly not too slow. I haven't verified the correctness of Ghosh's algorithm nor tested my own code very carefully, but it does correctly compute all the values listed on the OEIS A001055 page.

Here's the code; it uses my own little library of mathematical functions but anyone who is interested enough to try to run it at all is unlikely to have trouble finding or making implementations of the gcd and all_divisors functions (noting that all_divisors is only needed for _small_ numbers here, and trial division will be fine; this is also why I have implemented the ordinary integer partition function in the ridiculous way you can see below). My apologies in advance if the indentation gets messed up between me and the reader.

import functools, itertools
from fractions import Fraction as F
from g import gcd, all_divisors, factor_with_multiplicities

def intify(x):
    if x.denominator == 1: return x.numerator
    return x

def content(f):
    """GCD of all entries in f."""
    c = 0
    for a in f:
        c = gcd(c, a)
        if c == 1: return c
    return c

def lam(f):
    return intify(sum(F(1,d) for d in all_divisors(content(f))))

def smaller(f):
    return itertools.product(*(range(a+1) for a in f))

def sub(f,g):
    return tuple(a-b for (a,b) in zip(f,g))

def accumulate(a, coeff, f):
    for (j,b) in enumerate(f):
        a[j] += coeff*b

def p(n):
    """Number of partitions of n."""
    assert(n <= 20)
    return [1,1,2,3,5,7,11,15,22,30,42,56,77,101,135,176,231,297,385,490,627][n]
    # (ha ha)

@functools.cache
def mp(f):
    """Number of multiplicative partitions of a number whose prime factorization
    is product pj^f[j]; f is a tuple of nonnegative integers."""
    if sum(1 for a in f if a>0) <= 1: return p(sum(f))
    pf = [0 for a in f]
    for g in smaller(f):
        if all(a==0 for a in g): continue
        accumulate(pf, intify(lam(g)*mp(sub(f,g))), g)
        #print("  ", pf, g, lam(g), mp(sub(f,g)))
    pf = [intify(a) for a in pf]
    # pf should now be p(f) f
    result = sum(pf) // sum(f)
    #print(f, "->", result)
    return result

def mp1(n):
    return mp(tuple(r for (p,r) in factor_with_multiplicities(n)))

-- 
g

[David Radcliffe]

I get the same result in about 1 minute with the following code:

from functools import cache
from sympy import divisors
from math import isqrt

@cache
def a(n, m=2):
    return 1 + sum(a(n//d, d) for d in divisors(n) if m <= d <= isqrt(n))

I don't think that the 'easy' label is justified, unless someone finds a much more efficient algorithm.

David

To view this discussion visit https://groups.google.com/d/msgid/seqfan/93d4097b-83b8-4f4d-9917-39e1aedc8c94%40pobox.com.

[Gareth McCaughan]

On 22/05/2026 22:09, David Radcliffe wrote:
> I get the same result in about 1 minute with the following code:
>
> from functools import cache
> from sympy import divisors
> from math import isqrt
>
> @cache
> def a(n, m=2):
>     return 1 + sum(a(n//d, d) for d in divisors(n) if m <= d <= isqrt(n))

If that works that well (I expected it not to, evidently for no very
good reason, which is why I did something more complicated) then
probably there's a faster variation on the theme that only does one
factorization, operates entirely on the vector of exponents, and uses
something other than "compare the corresponding _numbers_" to do the
upper-bound thing with the second parameter.

--
g

[Christian Sievers]

This is essentially counting multiset partitions of multisets.

Here is my take using GP/PARI:

-----------------------------------------------

var(n) = eval(concat("'x",Str(n)))

s(p,k,n) = for(i=1,k,p=subst(p,var(i),var(i)^n)); p

a(n) = { my(f=factor(n)~[2,],
            o=prod(i=1,#f,1/(1-var(i))+O(var(i)^(f[i]+1)))-1,
            p=exp(sum(i=1,vecsum(f),s(o,#f,i)/i)));
         for(i=1,#f,p=polcoeff(p,f[i],var(i)));
         p
       }

-----------------------------------------------

It needs increased stack size, but then performs quite well:

? default(parisizemax, 80m)
  ***   Warning: new maximum stack size = 80003072 (76.297 Mbytes).
? #
   timer = 1 (on)
? a(13585283712000)
  *** subst: Warning: increasing stack size to 16000000.
  *** subst: Warning: increasing stack size to 32000000.
cpu time = 1,303 ms, real time = 1,303 ms.
%1 = 4258036618

Still not easy enough for "easy" if you ask me.

I'm sure it can still be improved, but I am already much longer awake than I wanted.

BTW, is there a better way to get an indexed variable than this var function?

All the best

Christian

[Christian Sievers]

[Had to reboot my computer.]

Sorry, the program works, but not in the indicated space and time, as it is not the version I used to get the output.

To get that, you have to replace "vecsum" with "vecmax" in the program.

Good night

[Antti Karttunen]

On Sat, May 23, 2026 at 7:09 AM Christian Sievers <g...@duvers.de> wrote:

Still not easy enough for "easy" if you ask me.

So, I'm still baffled, what is the OEIS definition of "easy"? Presumably, there should be a program/algorithm with which you can without too much pain compute a(n) for any somewhat arbitrarily large n, but in how short time? While "easy" should not be added to any sequence where you have to search for the next term, unless it is a "very easy search", that is guaranteed to yield results quite regularly, or? Note that A000040 has also keyword "easy". I guess the theory of computation doesn't help us much in setting any precise rules here?

Best,

Antti

[Neil Sloane]

I deleted the keyword "easy"!

Best regards

Neil 

Neil J. A. Sloane, Chairman, OEIS Foundation.

Also Visiting Scientist, Math. Dept., Rutgers University, 

Email: njas...@gmail.com

--
You received this message because you are subscribed to the Google Groups "SeqFan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to seqfan+un...@googlegroups.com.

To view this discussion visit https://groups.google.com/d/msgid/seqfan/CAB%2B0_%3D%3DnSBzUQO3CBpKvpj8BXQ08Y4sUEUf9snkby3V%2BTWR_aw%40mail.gmail.com.

[Allan Wechsler]

I was just thinking: probably nobody has thought hard about what the keyword "easy" means since Neil introduced the keyword. And then I thought "We should just ask Neil." 

I don't know if we do need a clear guideline. There seem to me to be two or three fairly fuzzy "easy zones", and some semantic wiggle-room having to do with what kind of sequence we are dealing with. There are several different, but related, tasks whose ease we might be thinking of when calling a sequence "easy".

• Given n, what is A(n)?
• Given k, is k in A?
• If k is in A, find n such that A(n) = k.
• Given the state of an A-producing algorithm that has just emitted A(n), how hard is it to keep running it to produce A(n+1)?

But I am merely woolgathering -- I'm not convinced we need clear guidelines about easiness. We usually know it when we see it, and that's probably good enough.

-- Allan

To view this discussion visit https://groups.google.com/d/msgid/seqfan/CAAOnSgT5bQ35kKXVJcmc8h4%2Bh4%2BS1%2BbMeTs%2B%3Dm_p9q9y3EJxfw%40mail.gmail.com.

[Antti Karttunen]

On Sat, May 23, 2026 at 7:16 PM Neil Sloane <njas...@gmail.com> wrote:

I deleted the keyword "easy"!

From A001055, but not from A000040. But surely A001055(13585283712000) is much easier to compute than A000040(13585283712000) ?!

Best,

Antti

To view this discussion visit https://groups.google.com/d/msgid/seqfan/CAAOnSgT5bQ35kKXVJcmc8h4%2Bh4%2BS1%2BbMeTs%2B%3Dm_p9q9y3EJxfw%40mail.gmail.com.

[Sean A. Irvine]

I tend to consider a sequence "easy" if we can quickly produce enough terms to do interesting follow on work. What I consider "quickly" or "enough" varies with the type of sequence, but say 10000 terms in under a second (by which A001055 does qualify as easy).  For "easy" it's more efficient to generate the terms than to download them from the b-file. I tend to think about it in terms of generating the sequence rather than a specific term of the sequence.

To view this discussion visit https://groups.google.com/d/msgid/seqfan/CADy-sGEt5dd6W-Qf1BBNyNOBQK05a76x0Sig7WQOXoXznRLukw%40mail.gmail.com.

[Charles Greathouse]

My writeup is here:

https://oeis.org/wiki/User:Charles_R_Greathouse_IV/Keywords/difficulty

To view this discussion visit https://groups.google.com/d/msgid/seqfan/CADy-sGEt5dd6W-Qf1BBNyNOBQK05a76x0Sig7WQOXoXznRLukw%40mail.gmail.com.

[Christian Sievers]

An obvious optimization for the substitution, lets our example run in less than 0.4s and without extra stack space.

Here is my latest version:

var(n) = eval(concat("'x",Str(n)))

s(p,f,n) = for(i=1,#f,my(v=var(i));p=subst(p+O(v^(1+f[i]\n)),v,v^n)); p

a(n) = { my(f=factor(n)[,2]~,
            o=prod(i=1,#f,1/(1-var(i))+O(var(i)^(f[i]+1)))-1,
            p=exp(sum(i=1,vecmax(f),s(o,f,i)/i)));
         for(i=1,#f,p=polcoef(p,f[i],var(i)));
         p
       }

The computed polynomial is symmetric under permutation of variables (except for the O(xi^k) terms) which

should allow to avoid repeated computations.

My remark about "easy" was just expressing my gut feeling, I didn't remember what I had read about this keyword before.

A40 being "easy" is a bit like me saying this problem is essentially multiset partition when my algorithm starts with a factorization,

as if that wasn't a hard problem!

Best,

Christian