Converting this recursive function definition into nested sums?

Here's the recursive function I'm using:

DH[n_, k_, s_] := Sum[Binomial[k, k - j] DH[n/m^j, k - j, m + 1], {m, s, n^(1/k)}, {j, 1, k}]
DH[n_, 0, s_] := 1


Now, on paper, if I expand out the recursion of DH by hand, I know that

DH[n,1,1] == n.


I also know, expanding by hand, that

DH[n,2,1] == 2 Sum[Floor[n/k], {k, 1, Floor[n^(1/2)]}] - Floor[n^(1/2)]^2


Can Mathematica do that for me, especially for larger values (like, say, DH[n,8,1] or DH[n,10,2])?

I would like to type in DH[n,2,1] or DH[n,3,1] and get those sums, but right now, I just get a recursion depth exceeded error and no answers.

The following is a note, not part of the original question

Computing DH[n,2,1] by hand, For Mr. Wizard:

Here's our starting formula

DH[n,2,1] = Sum[Binomial[k, k - j] DH[n/m^j, k - j, m + 1], {m, s, n^(1/k)}, {j, 1, k}]


apply s=1, k=2 and remove the {j,1,k} Sum by manually looping j from 1 to k, with k=2

DH[n,2,1] = Sum[ Binomial[2,1] DH[ n/m, 2-1,m+1 ] + Binomial[2,2] DH[n/m^2, 2-2, m+1],
{m,1,n^(1/2)}]


apply binomials, simplify trivial arithmetic

DH[n,2,1] = Sum[ 2 DH[n/m, 1,m+1] + DH[n/m^2,0,m+1], {m,1,n^(1/2)}]


Swap in DH[n,0,s] = 1

 DH[n,2,1] = Sum[ 2 DH[n/m, 1,m+1] + 1, {m,1,n^(1/2)}]


Inspection should make clear that DH[ n/m, 1, m+1 ] is Floor[n/m] - (m+1) + 1

DH[n,2,1] = Sum[ 2 (Floor[n/m]-m) + 1, {m,1,n^(1/2)}]


Distribute the x2

DH[n,2,1] = Sum[ 2 Floor[n/m] - 2 m + 1, {m,1,n^(1/2)}]


Extract the 1 from the sum

DH[n,2,1] = Floor[n^(1/2)] + Sum[ 2 Floor[n/m] - 2 m, {m,1,n^(1/2)}]


Use the triangle number formula (n)(n+1)/2 to pull -2 m out of the sum as -2 ( Floor[n^(1/2)] (Floor[n^(1/2)] + 1))/2

DH[n,2,1] = Floor[n^(1/2)] -2 ( Floor[n^(1/2)] (Floor[n^(1/2)] + 1))/2 + Sum[ 2Floor[n/m],
{m,1,n^(1/2)}]


Cancel the 2s outside the sum, and move the 2 inside the sum out front

DH[n,2,1] = Floor[n^(1/2)] - ( Floor[n^(1/2)] (Floor[n^(1/2)] + 1)) + 2 Sum[ Floor[n/m],
{m,1,n^(1/2)}]


Multiply the terms in parenthesis outside the sum

DH[n,2,1] = Floor[n^(1/2)] - Floor[n^(1/2)]^2 - Floor[n^(1/2)] + 2 Sum[ Floor[n/m],
{m,1,n^(1/2)}]


get rid of cancelling terms, and we are left with

DH[n,2,1] = -Floor[n^(1/2)]^2 + 2 Sum[ Floor[n/m], {m,1,n^(1/2)}]


...which is what I wanted to show. Check http://en.wikipedia.org/wiki/Divisor_summatory_function and you should see that is indeed the equations for D(n) when using the hyperbola method. That page mentions a generalized Divisor function, D_k(n). The k in that function corresponds to the k here - this recursive expression applies the hyperbola method to D_k(n) for any k, a whole number.

By way of comparison, I have simplified DH[x,4,2] to

1 - Floor[x^(1/4)]^4 +
Sum[
4 (Floor[x/(u^3)] + Floor[Floor[x/u]^(1/3)]^3)
- 6 Floor[ Floor[x/(u^2)]^(1/2)]^2 +
12 (Sum [Floor[x/(u^2 s)], {s, (u + 1),
Floor[Floor[x/(u^2)]^(1/2)]}] +
Sum[Floor[x/( u s^2)] - Floor[Floor[x/(u s)]^(1/2)]^2, {s, (u + 1),
Floor[Floor[x/u]^(1/3)]}]) +
24  Sum[ Floor[x/(u m s)], {s, u + 1, (x/u)^(1/3)}, {m, s + 1, (x/(u s))^(1/2)}],
{u, 2, x^(1/4)}]


I haven't gotten up the gumption to go higher, as it's a bit exhausting to work through - thus my question.

-
Humor me; would you include an example of "expanding by hand" DH[n,2,1]? –  Mr.Wizard Aug 4 '13 at 22:06
en.wikipedia.org/wiki/Divisor_summatory_function DH is a tidy recursive way of expressing the Hyperbola method for the Dirichlet Divisor problem (if s=1). DH[n,2,1] calculates D_2(n), from that link. –  Nathan McKenzie Aug 4 '13 at 22:10
That's really going the extra mile for my curiosity. Thank you! I wish I could up-vote a second time. –  Mr.Wizard Aug 4 '13 at 23:08

You can try either

DH[n_, k_, s_] :=
Sum[Sum[Binomial[k, k - j] DH[n/m^j, k - j, m + 1], {j, 1, k}], {m, s, n^(1/k)}]


or

DH[n_, k_, s_] :=
Sum[Sum[Binomial[k, k - j] DH[n/m^j, k - j, m + 1], {m, s, n^(1/k)}], {j, 1, k}]


Then you get

DH[n, 1, 1]
(* n *)


and

DH[n, 2, 1]
(* -n + 2 n HarmonicNumber[Sqrt[n]] *)


in both cases. However, the answer for DH[n, 2, 1] does not always agree with the answer you calculated by hand:

Table[-n + 2 n HarmonicNumber[Sqrt[n]] ==
2 Sum[Floor[n/k], {k, 1, Floor[n^(1/2)]}] - Floor[n^(1/2)]^2,
{n, 10}]
(* {True, False, False, True, False, False, False, False, False, False} *)


I'll leave you to figure out whether the error is in the code or the calculation.

However, while Mathematica does reduce the double sum to a single sum for larger values, it does not succeed in reducing the sum DH[n, 3, 1] to closed-form, nor DH[n, 8, 1] (DH[n, 10, 2] ran for several minutes and I aborted it). Perhaps someone else can help with that.

-
Thanks for the suggestion. Actually, if you put a Floor[] around n/m^j, the answer for DH[n,2,1] agrees with mine (except it's in a less simplified form), so that must've been the sticking point. It's still hanging on DH[n,3,1], though. But this gives me something to experiment with. –  Nathan McKenzie Aug 5 '13 at 0:13