# Can this be solved even faster?

So I would like to solve the following set of equation for $$m_i$$ given a set of $${M_m,N_m}$$.

$$m_1 +m_2 +m_3 +m_4 =M_m \\ |m_1| +|m_2| +|m_3| +|m_4| =N_m$$

All variables are integers. Also $$N_m \ge M_m$$ and their maximum value can reach up-to 30. I only need the total number of possible solution not the solutions themselves. So my first trivial attempt was to just use Solve

dimNM1[Nm_, Mm_] :=
Length[(Solve[m1 + m2 + m3 + m4 == Mm &&
Abs[m1] + Abs[m2] + Abs[m3] + Abs[m4] == Nm, {m1, m2, m3, m4}, Integers])]


My second slightly non-trivial attempt is the following:-

dimNM2[Nm_, Mm_] :=
Which[Nm === Mm,
Length[Partition[
Flatten[Permutations /@ IntegerPartitions[Nm, {4}, Range[0, Nm]]],
4]], True,
Module[{res},
res = Partition[
Flatten[Permutations /@ IntegerPartitions[Mm, {4}, Range[-Nm, Nm]]],
4];
Length[
Select[res, (Abs[#[[1]]] + Abs[#[[2]]] + Abs[#[[3]]] +
Abs[#[[4]]]) == Nm &]]]]


The second method is much faster than the first specially for $$N_m=M_m$$. But I would like to increase the speed further for $$N_m\ge M_m$$ case if possible.

dimNM1[2, 2] // AbsoluteTiming
(*{0.177768, 10}*)

dimNM2[2, 2] // AbsoluteTiming
(*{0.0000899056, 10}*)


So is there any other way to solve these equation faster?

• Note that N has built-in meanings. – Αλέξανδρος Ζεγγ Dec 9 '18 at 12:34
• OK I have changed it. – Hubble07 Dec 9 '18 at 12:46
• Nice problem. No need to generate candidates... see my reply. – ciao Dec 10 '18 at 8:13

ClearAll[num];

num[n_, m_] /; OddQ[n + m] = 0;
num[n_, n_] := Binomial[n + 3, 3];
num[n_, m_] /; OddQ[n] := With[{z = Ceiling[m/2]}, (5*n^2 + 3)/2 + 2 z - (2 z^2)];
num[n_, m_] /; EvenQ[n] := With[{z = Ceiling[m/2]}, (5*n^2 + 4)/2 - (2 z^2)];


Testing vs fastest answer here at writing (Henrik Schumacher):

stop = 100;

res = Table[{n, m, dimNM3[n, m]}, {n, 1, stop}, {m, 1, n}]; // AbsoluteTiming//First
res2 = Table[{n, m, num[n, m]}, {n, 1, stop}, {m, 1, n}]; // AbsoluteTiming//First

res == res2


169.203

0.0219434

True

Large cases are a non-issue:

num[123423456, 123412348] // AbsoluteTiming


{0.0000247977, 30468069908023290}

Some quick timings:

• Pretty impressive. Would you mind to elaborate where these formulas come from or at least to provide an (accessible) source? – Henrik Schumacher Dec 10 '18 at 8:49
• @HenrikSchumacher - I derived them, looking at a set of results: I recognized the pattern(s). Neat that the tetrahedral numbers and coordination sequences popped out. See e.g. Sloan, "Low-Dimensional Lattices VII: Coordination Sequences". – ciao Dec 10 '18 at 9:29
• Chapeaux for recognizing the patterns! =D – Henrik Schumacher Dec 10 '18 at 10:14
• @ciao - You Sir are a genius. Thank you. – Hubble07 Dec 10 '18 at 14:02
• Answers from ciao are generally great reads, +1. – Marius Ladegård Meyer Dec 11 '18 at 14:19

It is more efficient to first pick the integer partitions whose absolute values sum up to n before generating the permutations.

dimNM3[n_, m_] := Total[
Map[
Length@*Permutations,
Pick[#, Abs[#].ConstantArray[1, 4], n] &[
IntegerPartitions[m, {4}, Range[-n, n]
]
]
]
];

m = 20;
n = 40;
dimNM1[n, m] // AbsoluteTiming
dimNM2[n, m] // AbsoluteTiming
dimNM3[n, m] // AbsoluteTiming


{0.116977, 3802}

{0.995365, 3802}

{0.005579, 3802}

Sorry for not knowing much Mathematica, but I have a Python solution you might be able to follow. I'm putting this on the community wiki for anyone who wants to translate it.

def count_solutions(Nm, Mm):
firsthalves = dict()
for m1 in range(-Nm,Nm+1):
for m2 in range(-Nm,Nm+1):
m = m1+m2
n = abs(m1)+abs(m2)
key = (m,n)
if key in firsthalves:
firsthalves[key] += 1
else:
firsthalves[key] = 1

solutions = 0
for m3 in range(-Nm,Nm+1):
for m4 in range(-Nm,Nm+1):
m = m3+m4
n = abs(m3)+abs(m4)
key = (Mm-m, Nm-n)
if key in firsthalves:
solutions += firsthalves[key]
return solutions


This is a meet in the middle strategy. I enumerate all the possible $$m1,m2$$ combinations and record how many times each $$m1+m2,|m1|+|m2|$$ combination occurs in a dictionary.

Then I go through all the possible $$m3,m4$$ combinations and for each combination I calculate the necessary $$m1+m2,|m1|+|m2|$$ combination to make $$Mm,Nm$$, and I refer to the dictionary to find out how many $$m1,m2$$ combinations can make that.

The difference is that you go through the $$m1,m2$$ combination then the $$m3,m4$$ combinations, and the number of operations is roughly a square root of going through every $$m1,m2,m3,m4$$ combination. You should be able to solve for $$Nm = 1000,Mn = 0$$ in a few seconds.

A different approach, tied to @Hubble7's other question, that has the same speed as @ciao's answer. The key is in noting the sum of the negative numbers and sum of the nonnegative numbers are each fixed values, and so it is just a counting problem when we have 1 negative and 3 nonnegative terms, then 2 and 2, then 3 and 1. We can then use Mathematica's NumberOfCompositions[ ] function.

For values of Mn and Nn define

pos = (Nn+Mn)/2
neg = (Nn-Mn)/2


where posis the sum of the positive numbers in {m1, m2, m3, m4} and neg is the sum of the absolute value of the negatives ( so it is a positive number).

Now use Mathematica's NumberOfCompositions[n, k ] function which gives you the count of all of the ways to divide integer n into k terms, including 0 terms. If we want to find the number of compositions not including 0 terms we calculate NumberOfCompositions[n - k, k].

Note that for k=1, we have NumberOfCompositions[n, 1] = 1

If we have k negative terms, then we have Binomial[4,k] ways to arrange them. This is just {4, 6, 4} for k = {1, 2, 3}

So for values of neg and pos

perms = 4 * NumberOfCompositions[pos, 3]
+ 6  NumberOfCompositions[neg - 2, 2] NumberOfCompositions[pos, 2]
+ 4 NumberOfCompositions[neg - 3, 3]


And finally converting it into a function that accepts Mn and Nm ( while stealing some code from @ciao)

numNew[n_, m_] /; OddQ[n + m] = 0;
numNew[n_, n_] := Binomial[n + 3, 3];
numNew[n_, m_] := 4*NumberOfCompositions[(n + m)/2, 3]
+ 6*NumberOfCompositions[(n - m)/2 - 2, 2] NumberOfCompositions[(n + m)/2, 2]
+ 4*NumberOfCompositions[(n - m)/2 - 3, 3]


Check the timing against @ciao's answer above

num[123423456, 123412348] // AbsoluteTiming

{0.0000390021, 30468069908023290}

my function

numNew[123423456, 123412348] // AbsoluteTiming

{0.0000369493, 30468069908023290}

It is about as fast as @ciao's and also suggests (to me!) an approach to this question: Find position without iterating