Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Suppose I have a dimension formula (for a Lie algebra representation) that is as follows: $$ d(a,b) = {(a+1)(b+1)(a+b+2) \over 2} $$ Now consider the surface $F(a,b,n) = 0 = d(a,b) -n$ where $n \in \Bbb N$. Is there a Mathematica function that would allow me to find triples $(a,b,n) \in \Bbb N^3$ on this surface?

share|improve this question
up vote 12 down vote accepted

If what you really want is to find all natural numbers $a$ and $b$ on the surface for a given $n$, a search among the divisors of $2n$ will do it efficiently even for sizable values of $a$ and $b$. To do this, we find all divisors $k$ of $2n$ (Divisors), then all divisors of the quotient $2n/k$ (Divisors again), forming triples $(u,v,w)$ for which $uvw=2n$ (Flatten@Table). Then it remains only to screen for those of the form $u=a+b+2, v=b+1, w=a+1$ (Cases) and remove any duplicates (Union).

d[a_, b_] := (a + 1) (b + 1) (a + b + 2)/2;

 With[{n = d[119, 599]}, 
    Outer[{#1, #2, 2 n/(#1 #2)} &, First[k], Last[k]], 
     {k, {{#}, Divisors[2 n/#]} & /@ Divisors[2 n]}], 2], 
      {c_, b_, a_} /; c == a + b && b >= a :> {a - 1, b - 1}]]]]


($0.06$ seconds to find all four solutions for $n=d(119, 599) = 25\ 920\ 000$.)

Although the method works for any positive $n$, this example begins with an $n$ formed from a particular $a$ and $b$ in order to check that we at least get $(a,b)$ back in the output. (Try it on some factorials, such as $13!$ [no output] and $14!$ [one solution].)


To solve $n = (a+1)(b+1)(a+b+2)/2$ for a given $n$, write $x=a+b+2$, $y=b-a$ (which is non-negative with no loss of generality), whence we can recover $a = (x-y)/2-1$ and $b=(x+y)/2-1$, and note

$$y^2 = x^2 - 8n/x.$$

Therefore we can reduce a quadratic-time search to a linear-time search (in the number of divisors of $n$) by ranging over $x$, computing the right-hand side, and checking whether it is a square:

AbsoluteTiming[With[{n = 3240},
    Cases[Divisors[8 n], x_ /; x^3 >= 8 n 
      && Sqrt[x^2 - 8 n/x] \[Element] Integers :> {x, Sqrt[x^2 - 8 n/x]}]], 
    {x_, y_} /; EvenQ[x - y] :> {(x - y)/2 - 1, (x + y)/2 - 1}]]]


Because it is possible to get some spurious answers where $x$ and $y$ have opposite parity (for then $a$ and $b$ are only half-integral), those are screened out when converting back to $a$ and $b$.

The speedup gets better with larger $n$; e.g., the answers for $20!$ can be computed in $12$ seconds ($\{\{1270079,1422719\},\{484703,2935295\},\{355679,3525119\},\{235199,4432319\}\}$).

share|improve this answer
This is fantastic thanks very much whuber – Moderat May 23 '13 at 20:17
You're welcome. You shouldn't be so quick to accept this answer, though: this community does a great job of improving answers over the course of a day or two. Much better answers might appear soon :-). – whuber May 23 '13 at 20:41
I appreciate the explanation in this answer. +1 – RunnyKine May 23 '13 at 20:51
The equation is a cubic (parameterized by $n$), so it's possible a more efficient method to obtain all solutions can be found by computing generators of its group of integral points. But typically these groups are so small--I haven't found one larger than $12$ yet--as a practical matter it might not be worth the effort. – whuber May 23 '13 at 21:31

How about:

numPoints = 10;
dd[a_, b_, n_] := (a + 1) (b + 1) (a + b + 2)/2 - n
 dd[aa, bb, nn] == 0 && aa > 0 && bb > 0 && nn > 0, {aa, bb, 
  nn}, Integers, numPoints]

One could be slicker, but this works.

share|improve this answer

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.