# How to speed up integers finding function?

Playing with integers I looked for a way of finding, let's say, $$4$$-digits positive integers such that $$a,b,a+b$$ had the same digits, like $$1089 + 8019 = 9108$$ I am a newbie and I used this function

sd[a_, b_] :=
If[Mod[a, 9] != 0 || Mod[b, 9] != 0 , False,
Sort[IntegerDigits[a + b]] == Sort[IntegerDigits[a]] &&
Sort[IntegerDigits[a]] == Sort[IntegerDigits[b]]]


Then I used the function in this way

Select[Flatten[
Table[{h, k, sd[h, k]}, {h, 1000, 10000}, {k, h, 10000}], 1], #[[3]] &]


But it took ages to give the output.

Is there a way to speed up this procedure?

• Shouldn't the iterators in the table be Table[..., {h, 1000, 10000 - 1}, {k, h, 10000 - h - 1}]? Oct 16, 2020 at 14:27
• The fastest way to find integers is to not lose them in the first place. Just saying. Oct 17, 2020 at 16:32

ClearAll[pairS]

pairS[n_] := SortBy[First] @
Apply[Join] @
KeyValueMap[Function[{k, v},
Select[k == Sort@IntegerDigits@Total@# &]@Subsets[v, {2}]]] @
GroupBy[Sort@*IntegerDigits] @
(999 + 9 Range[10^(n - 1)])


Examples:

 pairS[4] // AbsoluteTiming // First

0.0445052

pairS[5] // AbsoluteTiming // First

1.19877

Multicolumn[pairS[4], 5]


Length @ pairS[5]

673

pairS[5] // Short[#, 7] &


An aside: A slower graph-based method: get the edge list of a graph where the numbers $$a$$ and $$b$$ are connected if $$a$$, $$b$$ and $$a+b$$ have the same integer digits.

relation = Sort[IntegerDigits @ #] == Sort[IntegerDigits @ #2] ==
Sort[IntegerDigits[# + #2]] &;

relationgraph = RelationGraph[relation, 999 + 9 Range[10^(4 - 1)]];

edges = EdgeList @ relationgraph;

List @@@ edges == pairS[4]

True

Subgraph[relationgraph, VertexList[edges],
GraphLayout -> "MultipartiteEmbedding",
GraphStyle -> "VintageDiagram", ImageSize -> Large]


• I'll spend next week studying your answer. Thank you Oct 16, 2020 at 15:02
• On my laptop the $5$-digits version took almost 6 minutes. I wrote a simple program in Pascal and it took about 28 seconds. Can you explain in simple words the reason why this happens? Oct 16, 2020 at 15:10
• @Raffaele I think the main problem is generating all the possible candidates beforehand, which gets very memory-intensive. Once the RAM on your computer fills up and it becomes necessary to swap memory to disk, things slow down dramatically. Oct 16, 2020 at 15:15
• I'll spend next TWO weeks studying your code. TY Oct 16, 2020 at 18:27
• Maybe you should correct 999 + 9 Range[10^(n - 1)] in 10^(n - 1) + 9 Range[10^(n - 1)] Oct 16, 2020 at 19:23

https://oeis.org/A331468

The numbers of 3-digit to 8-digit triples are: 1, 25, 648, 17338, 495014, and 17565942.

Approach 1, more concise

Clear[search];
search[n_] :=
Join @@ Table[With[{s = Subsets[a, {2}]},
Pick[s, Boole@MemberQ[a, Total@#] & /@ s, 1]],
{a, GatherBy[Select[Range[10^(n - 1), 10^n - 1], Divisible[#, 9] &],
Sort@*IntegerDigits]}];

search[4] // Length // AbsoluteTiming
search[5] // Length // AbsoluteTiming
search[6] // Length // AbsoluteTiming


{0.0210189, 25}
{0.212638, 648}
{9.23615, 17338}

Approach 2, more efficient

Clear[cf]
cf = Compile[{{n, _Integer}, {A, _Integer, 2}},
Module[{nums, ni, nj, B = InternalBag[Most@{0}]},
Do[
nums = Permutations[a]. 10^Range[n - 1, 0, -1];
Do[
ni = nums[[i]];
nj = nums[[j]];
If[ni + nj > 10^n || ni < 10^(n - 1), Break[]];
Do[If[ni + nj == k, InternalStuffBag[B, {ni, nj, k}, 1]; Break[]]
, {k, nums}]
, {i, Length@nums}, {j, i + 1, Length@nums}]
, {a, A}];
InternalBagPart[B, All]
], CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

n = 4;
AbsoluteTiming[
digits = Select[# - Range[n] & /@ Subsets[Range[9 + n], {n}], Divisible[Total@#, 9] &];
Length[ans = Partition[cf[n, digits], 3]]
]


For n=4,5,6

{0.0014472, 25}
{0.0094707, 648}
{0.802517, 17338}

ClearAll[pairS]
pairS[n_] :=
Apply[Join]@ KeyValueMap[Function[{k, v},
Select[k == Sort@IntegerDigits@Total@# &]@Subsets[v, {2}]]]@
GroupBy[Sort@*IntegerDigits]@(10^(n - 1) - 1 + 9 Range[10^(n - 1)])

pairS[4] // Length // AbsoluteTiming
pairS[5] // Length // AbsoluteTiming
pairS[6] // Length // AbsoluteTiming


{0.0362128, 25}
{0.945485, 648}
{40.879, 17338}

Approach 3, Using LibraryLink, 10 times faster again

src="#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <WolframLibrary.h>

void solve(int *a, int level, int n, int start, std::vector<int> &res) {
if (level == n) {
if (std::accumulate(a, a + n, 0) % 9 == 0) {
std::vector<int> vec;
do {
if (a[0] != 0)
vec.push_back(std::accumulate(a, a + n, 0,
[&](int x, int y) { return 10 * x + y; }));
} while (std::next_permutation(a, a + n));
int &max = *std::max_element(vec.begin(), vec.end());
for (int i = 0; i < vec.size(); ++i) {
for (int j = i + 1; j < vec.size(); ++j) {
int ni = vec[i], nj = vec[j];
if (ni + nj > max)
break;
if (std::binary_search(vec.begin(), vec.end(), ni + nj)) {
res.push_back(ni);
res.push_back(nj);
res.push_back(ni+nj);
}
}
}
}
return;
}
for (int i = start; i <= 9; i++) {
a[level] = i;
solve(a, level + 1, n, i, res);
}
}

EXTERN_C DLLEXPORT int func(WolframLibraryData libData, mint Argc, MArgument *Args,
MArgument Res) {
MTensor out;
mint *out_data;
int n = MArgument_getInteger(Args[0]);
std::vector<int> res;
int *a = new int[n];
solve(a, 0, n, 0, res);
mint len = res.size();
int err = libData->MTensor_new(MType_Integer, 1, &len, &out);
out_data = libData->MTensor_getIntegerData(out);
std::copy(res.begin(), res.end(), out_data);
MArgument_setMTensor(Res, out);
return 0;
}
";

Needs["CCompilerDriver"];
\$CCompiler={"Compiler"->CCompilerDriverGenericCCompilerGenericCCompiler,
"CompilerInstallation"->"C:/msys64/mingw64","CompilerName"->"g++.exe",
"CompileOptions"->"-O3"};

Needs["CCompilerDriver"];
lib=CreateLibrary[src,"func"];

Table[n -> AbsoluteTiming[Partition[func[n], 3] // Length], {n, 2, 8}] // Column


• Thanks for posting this! I was wondering about (considering investigating) the use of Compile for this problem. Oct 17, 2020 at 17:12
• @AntonAntonov You are welcome. Oct 17, 2020 at 17:37
• If I'm correct, The numbers of 3-digit to 10-digit triples are: 1, 25, 648, 17338, 495014, 17565942, 717564880, 30694477548. How is the growth like?
– l4m2
May 22, 2023 at 7:53

But it took ages to give the output.

It took ~170 seconds on my computer; with ParallelTable it took ~97 seconds.

I assume two-times speed-up is not good enough, but it was very easy to get it.

Divide the numbers from 1000 to 9999 into a few hundred sets of integers that have the same digits, for example [1234, 1243, 1324, 1342, 1423, 1432 ... ]. Then a and b must be in the same set, and a+b must be in that set as well. So you loop over the 400 or so sets S of integers, then iterate over all elements a < 5000 of the set S, iterate b over all elements of the set S with a ≤ b ≤ 9999-a, and then check if a+b is an element of S as well. Should take milliseconds.

Maybe out of slope...

Since this range is kind of huge. So use Python's Api maybe a better choice?

ExternalEvaluate["Python", "[(i, j, i+j)for i in range(1000, 9999) for j in range(i, 9999-i)
if sorted(str(i)) == sorted(str(j)) == sorted(str(i+j))]"] // AbsoluteTiming

{27.2873, {{1089, 8019, 9108}, {1089, 8091, 9180}, {1269, 1692,
2961}, {1467, 6147, 7614}, {1467, 6174, 7641}, {1476, 4671,
6147}, {1503, 3510, 5013}, {1530, 3501, 5031}, {1746, 4671,
6417}, {2385, 2853, 5238}, {2439, 2493, 4932}, {2502, 2520,
5022}, {2538, 3285, 5823}, {2691, 6921, 9612}, {2853, 5382,
8235}, {3285, 5238, 8523}, {4095, 4950, 9045}, {4095, 5409,
9504}, {4392, 4932, 9324}, {4590, 4950, 9540}, {4599, 4995,
9594}, {4698, 4986, 9684}, {4797, 4977, 9774}, {4896, 4968,
9864}, {4959, 4995, 9954}}}
`

costs 27s

Contrast to origin code which takes 233.128s on my PC.

• (+1) but I was seriously contemplating to give this answer (-1). :) Oct 16, 2020 at 17:29
• This is Python, not Mathematica Oct 16, 2020 at 21:00
• @Raffaele "This is Python, not Mathematica" -- yes, but two things have to be said: 1. This solution uses Python called within Mathematica with native (to Mathematica) input and output. 2. Your (original) post finished with the question "Is there a way to speed up this procedure?" and this answer satisfies it. Oct 17, 2020 at 0:05