Hi,
This has the recurrence relation
T(n,k) = T(n-1,k-1) + T(n-k,k) OR
T(n,k) = sum(j=1,k,T(n-k,j))
With some initial/boundary conditions:
=1 when (k==1) || (k==n) || ((k==0) && (n==0))
=0 when (k>n) || (n != 0 && k < 1)
So I have some PARI code already which can do this(and indeed there is some listed in the link, which I started with), but I'm trying my best to make it more efficient for very large n values, without necessarily memoizing the entire table. One thing that can help speed this up is to have polynomial equations for specific small k values.
I have adapted some equations for some small k values which are:
T(n,1)=1
T(n,2)=floor(n/2))
T(n,3)=floor((n^2+6)/12)
T(n,4)=floor((n^3+3*n^2-9*n*(n%2)+72)/144)
T(n,5)=floor((n^4+10*(n^3+n^2)-75*n-45*n*(-1)^n+1440)/2880)
I actually don't know how these equations are determined beyond k=2, but I found them from various OEIS pages. My ultimate goal is to implement such a function in C using GMP or FLINT arbitrary precision libraries, but I'm using PARI as testing ground for getting the formulas right.
I use the floor function and add half the denominator for these(as opposed to calling round(...), which is how they were originally written), since this will translate better into integer division when I write the code for GMP. I'm also slightly skeptical of the rounding in these formula and I'm not certain that they will hold true (with error < 1) for any arbitrarily large n. Does the math always work out correctly that way?
The oeis page provides some generating functions for the sequence, but I just can't grasp how to make use of these:
G.f. for k-th column: x^k/(Product_{j=1..k} (1-x^j))
G.f.: A(x, y) = Product_{n>=1} 1/(1-x^n)^(P_n(y)/n), where P_n(y) = Sum_{d|n} eulerphi(n/d)*y^d
G.f.: G(t,x) = -1 + 1/Product_{j>=1} (1-t*x^j)
G.f.: -1 + e^(F(x,z)), where F(x,z) = Sum_{n >= 1} (x*z)^n/(n*(1 - z^n)) is a g.f. for A126988
Is there a way I can use PARI to help generate the formulas for specific values of k as above, by using these G.f.'s?
I've also been working on attacking the "other side" of this problem, using shortcuts for specific cases where ak+b>=n
Given the functions:
p(n) = numbpart(n) // total partitions of n
and
A000070(n)= sum(k=1,n,p(k))
And the equivalences:
p(n) == sum(k=1,n,T(n,k))
T(n,k) == sum(j=1,k,T(n-k,j))
It follows that:
T(n,k) == p(n-k) - sum(j=k+1,n, T(n-k,j))
I call the sum which is subtracted from p(n-k) the "anti-children" of T(n,k)
This can be applied to multiple times to get a formula which doesn't require recursively evaluating T(n,k).
Each successive step for ak+b>=n is defined such that the anti-children are guaranteed to satisfy the previous inequality:
k > n: T(n,k) = 0
2k >= n: T(n,k) = p(n-k)
3k+2>=n: T(n,k) = p(n-k) - A000070(n-2*k-1)
4k+5>=n: T(n,k) = p(n-k) - A000070(n-2*k-1)+Sum_{j=k+1..floor((n-k-1)/2)} A000070(n-k-2*j-1))
4k+5< n: T(n,k) = p(n-k) - A000070(n-k-floor((n-k)/3)) + (Sum_{j=floor((n-k)/3)..floor((n-k-1)/2)} A000070(n-k-2*j-1) ) - Sum_{j=k+1..floor((n-k)/3)-1} T(n-k,j)
These cases were already on the oeis page up to 3k+2>=n, then I spent quite some time working out the formula for 4k+5>=n by hand, (and verifying against a printed table, and fixing mistakes). I might try to continue this process (next would be 5k+9>=n, 6k+14>=n, 7k+20>=n, etc.) but it gets increasingly difficult for me.
Alternatively, if anyone understands how Hardy-Ramanujan-Rademacher equation works for p(n), could please tell me if its possible to apply this somehow to the restricted "n into k parts" problem, that would be amazing. The FLINT library uses this HRR for partitions of n, and is apparently incredibly efficient, but doesn't provide a function for "n into k" parts
Or if you have thoughts about any other ways to approach this problem, I would be interested to hear.
Thank you,
Hans Loeblich