12
$\begingroup$

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?

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

enter image description here

$\endgroup$
  • 3
    $\begingroup$ Pretty impressive. Would you mind to elaborate where these formulas come from or at least to provide an (accessible) source? $\endgroup$ – Henrik Schumacher Dec 10 '18 at 8:49
  • 3
    $\begingroup$ @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". $\endgroup$ – ciao Dec 10 '18 at 9:29
  • 2
    $\begingroup$ Chapeaux for recognizing the patterns! =D $\endgroup$ – Henrik Schumacher Dec 10 '18 at 10:14
  • 3
    $\begingroup$ @ciao - You Sir are a genius. Thank you. $\endgroup$ – Hubble07 Dec 10 '18 at 14:02
  • 1
    $\begingroup$ Answers from ciao are generally great reads, +1. $\endgroup$ – Marius Ladegård Meyer Dec 11 '18 at 14:19
13
$\begingroup$

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}

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

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