7
$\begingroup$

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?

Thanks in advance

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

5 Answers 5

12
$\begingroup$
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]

enter image description here

Length @ pairS[5] 
673
pairS[5] // Short[#, 7] &

enter image description here

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]

enter image description here

$\endgroup$
8
  • 2
    $\begingroup$ I'll spend next week studying your answer. Thank you $\endgroup$
    – Raffaele
    Oct 16, 2020 at 15:02
  • $\begingroup$ 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? $\endgroup$
    – Raffaele
    Oct 16, 2020 at 15:10
  • 1
    $\begingroup$ @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. $\endgroup$ Oct 16, 2020 at 15:15
  • 2
    $\begingroup$ I'll spend next TWO weeks studying your code. TY $\endgroup$
    – Raffaele
    Oct 16, 2020 at 18:27
  • 1
    $\begingroup$ Maybe you should correct 999 + 9 Range[10^(n - 1)] in 10^(n - 1) + 9 Range[10^(n - 1)] $\endgroup$
    – Raffaele
    Oct 16, 2020 at 19:23
5
$\begingroup$

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 = Internal`Bag[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, Internal`StuffBag[B, {ni, nj, k}, 1]; Break[]]
       , {k, nums}]
      , {i, Length@nums}, {j, i + 1, Length@nums}]
     , {a, A}];
    Internal`BagPart[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}

Compare with kglr's answer

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"->CCompilerDriver`GenericCCompiler`GenericCCompiler,
"CompilerInstallation"->"C:/msys64/mingw64","CompilerName"->"g++.exe",
"CompileOptions"->"-O3"};

Needs["CCompilerDriver`"];
LibraryUnload["func"];
lib=CreateLibrary[src,"func"];
func=LibraryFunctionLoad[lib,"func",{Integer},{Integer,1}];

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

enter image description here

$\endgroup$
3
  • 1
    $\begingroup$ Thanks for posting this! I was wondering about (considering investigating) the use of Compile for this problem. $\endgroup$ Oct 17, 2020 at 17:12
  • 1
    $\begingroup$ @AntonAntonov You are welcome. $\endgroup$
    – chyanog
    Oct 17, 2020 at 17:37
  • $\begingroup$ 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? $\endgroup$
    – l4m2
    May 22, 2023 at 7:53
2
$\begingroup$

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.

enter image description here

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ (+1) but I was seriously contemplating to give this answer (-1). :) $\endgroup$ Oct 16, 2020 at 17:29
  • 2
    $\begingroup$ This is Python, not Mathematica $\endgroup$
    – Raffaele
    Oct 16, 2020 at 21:00
  • 1
    $\begingroup$ @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. $\endgroup$ Oct 17, 2020 at 0:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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