Tag List is Protected Error in simple script

Question

I want to have a simple script that gives me all numbers $k \in \mathbb{N}$ in some range such that product of all their divisors is equal to $k^3$. Here is the code that I wrote:

  findDivisors[n_] := Module[
   {i, answer},
   answer = {};
   For[i = 1, i <= n, i++,
    If[Apply[Times, Divisors[i]] == i^3, Append[answer, i]]];
   Return[answer];
   ];

So it must check all integers from $1$ to $n$ and gives me the list "answer" which elements satisfy the desired property. But when I run it the output is

SetDelayed::write: Tag List in {}[n_] is Protected.

So where is my mistake?

Answer 1

I get no such error from your code but I also get no useful output because you used Append where you needed AppendTo. See Assuming commands will have side effects when they don't.

Please try:

ClearAll[findDivisors]

findDivisors[n_] :=
  Module[{i, answer},
   answer = {};
   For[i = 1, i <= n, i++,
    If[Apply[Times, Divisors[i]] == i^3, AppendTo[answer, i]]];
   answer
  ]

findDivisors[100]

{1, 12, 18, 20, 28, 32, 44, 45, 50, 52, 63, 68, 75, 76, 92, 98, 99}

A more idiomatic way to write that is:

findDivisors2[n_] :=
  Table[
    If[Times @@ Divisors[i] == i^3, i, ## &[]],
    {i, n}
  ]

Or taking additional liberties with the translation:

findDivisors3[n_] := Cases[Range @ n, i_ /; Times @@ Divisors[i] == i^3]

Reference:

• How to avoid returning a Null if there is no "else" condition in an If contruct
• Alternatives to procedural loops and iterating over lists in Mathematica

---

Performance

Using LLlAMnYP's observation that each output number besides one has exactly six divisors here is a function targeting performance making use of DivisorSigma:

findDivisors4[n_] :=
  Pick[#, DivisorSigma[0, #], 6] & @ Range @ n // Prepend[1]

Needs["GeneralUtilities`"]

BenchmarkPlot[{findDivisors, findDivisors2, findDivisors3, findDivisors4}, # &]

Though this does not improve on the computational order of the other methods is it constantly three times faster.

Answer 2

Like I commented, every single number required by the OP has exactly six divisors (including one and itself, with the trivial exception of 1). Consequentially, there are only two possible options: either the number is a fifth power of a prime: $p^5$ and has divisors $p^{0..5}$ whose product is, of course, $p^{15} = {(p^5)}^3$; otherwise it is the product $p_i^2 p_j$ and has divisors $1, p_i, p_j, p_i p_j, p_i^2, p_i^2 p_j$. Their product is of course $p_i^6 p_j^3 = {(p_i^2 p_j)}^3$.

I now want to generate all possible tuples of the form {i, i, j} for i<=n and j<=n (and i != j). These kind of tasks are best met by built-in functions like Subsets, but since the first number is repeated, a customized compiled function to suit my needs should also work very well.

Code dump:

getTuples = 
 Compile[{{n, _Integer}}, 
  Block[{out = Table[0, {n * (n - 1)}, {3}], i = 1, j = 1, k = 1},
   For[i = 1, i <= n, i++,
    For[j = 1, j < i, j++; k++,
     out[[k, {1, 2}]] = Table[i, {2}]; out[[k, 3]] = j];
    For[j = i + 1, j <= n, j++; k++,
     out[[k, {1, 2}]] = Table[i, {2}]; out[[k, 3]] = j]
    ]; out]]

to actually get the required numbers we run the following:

getTuples[10] // Prime // Times @@@ # &

We lack here the fifth powers of primes though. To remedy, I define a little helper function:

getNumbers[n_] := 
  Sort[{1}~Join~((Prime[Range[n]])^5)~
    Join~(getTuples[n] // Prime // Times @@@ # &)]

It gets me 1000001 of such numbers in 3 seconds, however there is no guarantee that they are the first 1000001; since $p_i^2 p_j$ for $j = 1, i = 1000$ is present, but the likely much smaller case of, say, $j = 1001, i = 2$ is not.

Obligatory sample output:

getNumbers[10]
(* {1, 12, 18, 20, 28, 32, 44, 45, 50, 52, 63, 68, 75, 76, 92, 98, 99, \
116, 117, 147, 153, 171, 175, 207, 242, 243, 245, 261, 275, 325, 338, \
363, 425, 475, 507, 539, 575, 578, 605, 637, 722, 725, 833, 845, 847, \
867, 931, 1058, 1083, 1127, 1183, 1421, 1445, 1573, 1587, 1682, 1805, \
1859, 2023, 2057, 2299, 2523, 2527, 2645, 2783, 2873, 3125, 3179, \
3211, 3509, 3703, 3757, 3887, 3971, 4205, 4693, 4901, 5491, 5819, \
5887, 6137, 6647, 6877, 8303, 8381, 8993, 9251, 10051, 10469, 10933, \
14297, 15341, 15979, 16807, 19343, 161051, 371293, 1419857, 2476099, \
6436343, 20511149} *)

Observe how it fails though:

ListPlot[{getNumbers[10], getNumbers[100][[;; 101]]}, PlotRange -> {Full, {0, 1000}}]

The output of the first 19 primes matches for both lists, but then getNumbers[10] shoots up, missing several suitable values. There may be ways to improve this.