Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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], 

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], 

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], 

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], 

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 +
      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.

share|improve this 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

1 Answer 1

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)}]


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 *)


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.

share|improve this answer
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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.